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Abstract 

Transformation  optics  has  shown  the  ability  to  cloak  an  object  from  incident 
electromagnetic  radiation  is  theoretically  possible.  However,  the  constitutive  param¬ 
eters  dictated  by  the  theory  are  inhomogeneous,  anisotropic,  and,  in  some  instances, 
singular  at  various  locations.  In  order  for  a  cloak  to  be  practically  realized,  sim¬ 
plified  parameter  sets  are  required.  However,  the  simplified  parameters  result  in  a 
degradation  in  the  cloaking  function. 

Constitutive  parameters  for  simplified  two-dimensional  cylindrical  cloaks  have 
been  developed  with  two  specific  material  property  constraints.  It  was  initially  be¬ 
lieved  satisfying  these  two  constraints  would  result  in  the  simplified  cylindrical  cloaks 
satisfying  the  same  wave  equation  as  an  ideal  cloak.  Because  of  this  error,  the  sim¬ 
plified  two-dimensional  cylindrical  cloaks  were  not  perfect.  The  error  in  the  initial 
derivation  of  the  original  simplified  parameter  sets  was  noted  in  the  published  litera¬ 
ture.  However,  no  analysis  was  done  to  determine  all  material  parameter  constraints 
to  enable  a  perfect  two-dimensional  cylindrical  cloak.  This  research  developed  a 
third  constraint  on  the  material  parameters.  It  was  shown  as  the  material  param¬ 
eters  better  satisfy  this  new  equation,  a  two-dimensional  cylindrical  cloak’s  hidden 
region  is  better  shielded  from  incident  radiation.  Additionally,  a  novel  way  to  derive 
simplified  material  parameters  for  two-dimensional  cylindrical  cloaks  was  developed. 
A  Taylor  series  expansion  dictated  by  the  new  constraint  equation  led  to  simplified 
cloaks  with  significantly  improved  scattering  width  performances  when  compared  to 
previous  published  results. 

During  the  course  of  this  research,  it  was  noted  all  cloak  simulations  are  per¬ 
formed  using  finite  element  method  (FEM)  based  numerical  methods.  While  accu¬ 
rate,  FEM  methods  can  be  computationally  intensive  and  time  consuming.  A  Green’s 
function  was  used  to  accurately  calculate  scattering  widths  from  a  two-dimensional 
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cylindrical  cloak  with  a  perfect  electrically  conducting  inner  shell.  Significant  time 
improvements  were  achieved  using  the  Green’s  function  compared  to  an  FEM  solution 
particularly  as  the  computational  domain  size  is  increased. 

Finally,  cloaks  are  physically  realized  using  metamaterials.  Design  of  metama¬ 
terials  has  typically  been  done  empirically.  Shifts  in  S-parameter  measurements  and 
the  resulting  extracted  constitutive  parameters  are  used  to  determine  the  impact  to 
resonant  regions  due  to  various  geometries.  A  new  way  to  design  and  possibly  opti¬ 
mize  unit  cell  metamaterials  was  investigated  using  an  eigendecomposition  method 
to  identify  unit  cell  resonances.  Different  structures  were  shown  to  have  different 
resonances,  and  control  of  the  resonant  locations  can  lead  to  optimum  designs. 
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Electromagnetic  Field  Control 


and 

Optimization  Using  Metamaterials 

I.  Introduction 

Radar  cross  section  (RCS)  reduction  has  been  a  goal  of  scientists  and  engineers 
since  the  first  major  uses  of  radar  in  World  War  II.  A  wide  body  of  knowledge  exists  on 
passive  techniques  used  to  control  scattering  due  to  incident  electromagnetic  energy. 
These  techniques  can  be  divided  into  two  main  sub-categories:  shaping  and  radar  ab¬ 
sorbing  materials  (RAM).  The  effectiveness  of  both  categories  are  typically  dependent 
on  the  frequency,  incident  angle,  and  polarization  of  the  illuminating  energy. 

The  goal  of  shaping  is  to  scatter  the  incident  energy  from  the  target  such  that  the 
amount  of  energy  returned  toward  the  radar  is  minimized.  This  type  of  RCS  control 
has  proven  to  be  very  effective  for  monostatic  radars  where  the  transmit  and  receive 
antennas  are  collocated.  For  a  threat  aircraft,  the  attack  profile  can  be  controlled 
such  that  a  small  range  of  target  angles  will  be  presented  to  the  radar.  Shaping  can 
be  used  to  reduce  the  RCS  at  these  angles  and  increase  the  aircraft’s  stealthiness. 
However,  a  rule  of  thumb  for  shaping  is  a  reduction  in  the  RCS  at  one  aspect  angle  is 
always  accompanied  by  an  increase  at  another  [49] .  Consider  a  two-dimensional  RCS 
(i.e.  echowidth).  If  all  360°  of  measurement  angles  are  equally  important,  shaping  will 
reduce  the  echowidth  at  one  angle  while  increasing  it  at  another  angle  (or  several  other 
angles).  In  some  instances,  this  is  acceptable.  As  an  example,  the  technique  of  lobe 
width  control  allows  the  RCS  to  balloon  in  certain  sectors  where  significant  scattered 
energy  does  not  impact  the  desired  result.  In  other  instances,  lobe  width  control  is 
an  unacceptable  way  to  control  RCS.  It  is  easy  to  see  a  low-RCS  threat  designed  to 
act  against  a  monostatic  radar  will  have  a  significantly  reduced  capability  against  a 
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bistatic  radar,  whether  the  bistatic  radar  is  an  integrated  network  of  transmitters  and 
receivers  or  if  the  radar  is  simply  making  use  of  transmitters  of  opportunity. 

The  goal  of  RAM  is  to  convert  the  incident  electromagnetic  energy  into  heat, 
thereby  reducing  the  amount  of  scattered  electromagnetic  energy  capable  of  being 
detected  by  the  radar  sensors.  Note  RAM  cannot  reduce  the  forward  scatter  of  an 
object,  but  it  is  very  effective  in  controlling  the  back  scatter.  This  can  be  done  with 
either  an  electric  or  magnetic  loss  tangent  in  the  RAM  material.  Narrow  band  RAM 
coatings,  such  as  the  Salisbury  screen  and  Dallenbach  layer,  have  been  used  since 
the  1950’s  [94],  Modern  radar  systems  span  a  wide  range  across  the  electromagnetic 
spectrum,  with  most  operating  between  220  MHz  -  35  GHz;  however,  over-the-horizon 
and  millimeter  wave  radars  operate  outside  this  range  [87].  Many  radars  also  have 
a  wide  operational  bandwidth.  Therefore,  wide  band  RAM  is  very  desirable.  The 
first  broadband  absorber  was  a  Jaumann  absorber,  which  can  be  thought  of  as  a 
multilayered  Salisbury  screen  [94],  Another  type  of  broadband  RAM  are  the  carbon- 
loaded  foam  absorbers  used  in  anechoic  chambers  to  limit  the  scattered  energy  from 
surfaces  other  than  the  device-under-test  (DUT)  [48].  Typical  RAM  employed  on 
modern  aircraft  is  some  type  of  iron  ball  paint.  The  paint  contains  tiny  spheres 
coated  with  carbonyl  iron  or  ferrite.  Incident  electromagnetic  energy  interacts  with 
these  spheres,  resulting  in  the  electromagnetic  energy  being  converted  to  heat  [3]. 

There  are  significant  implications  when  using  RAM.  First,  most  are  toxic  to 
some  degree.  During  the  first  Gulf  War,  maintenance  crews  noted  a  large  number  of 
dead  bats  in  the  hangars  where  the  F-117  was  kept.  Their  deaths  were  attributed  to 
long  exposure  to  the  RAM  coupled  with  a  lack  of  ventilation  [1] .  Additionally,  RAM 
coatings  require  precise  application  methods,  as  the  coating  thickness  and  smoothness 
must  be  uniform  across  the  surface  of  the  substrate.  The  application  process  typi¬ 
cally  involves  robotic  sprayers  that  can  accurately  control  the  coating  thickness  [3]. 
Furthermore,  the  applied  coatings  require  strict  constitutive  parameter  tolerances  as 
well  as  uniformity  in  order  to  achieve  the  desired  result.  Therefore,  costs  increase 
drastically  when  working  with  RAM.  Also,  any  type  of  RAM  coating  increases  an 
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object’s  weight.  For  aircraft,  weight  increases  can  have  significant  impact  on  perfor¬ 
mance.  RAM  is  not  simply  a  covering  that  can  be  easily  applied  to  an  aircraft  or 
other  body  to  reduce  its  RCS.  Rather,  working  with  RAM,  from  its  manufacture  to 
application,  is  a  technically  detailed,  costly  process. 

When  used  together,  shaping  and  RAM  make  effective  RCS  reduction  tools,  but 
the  limits  of  their  effectiveness  are  being  approached  in  the  field  of  RCS  measurements. 
These  limits  are  particularly  noticeable  when  the  DUT  has  a  low  RCS.  The  low  RCS 
makes  it  difficult  to  differentiate  radar  returns  resulting  from  the  DUT  and  other 
objects  inherent  in  the  measurement  system.  This  limitation  is  discussed  in  the 
following  section. 

1.1  Long  Term  Problem  Statement 

RCS  measurements  are  obtained  using  test  ranges.  Static  test  ranges  can  be 
either  indoor  or  outdoor,  with  each  having  its  positives  and  negatives.  For  physically 
large  targets,  an  outdoor  test  range  is  typically  required.  The  ideal  outdoor  test 
range  has  minimal  background  signals  and  little  to  no  secondary  scattering  sources. 
This  allows  the  measured  RCS  to  be  as  close  to  the  actual  RCS  as  possible.  To 
avoid  ground  bounce  interaction,  the  target  is  mounted  a  considerable  distance  off 
the  ground.  A  metal  pylon  is  often  used  as  the  primary  support  structure.  Other 
support  systems  exist,  such  as  foam  columns  and  string  support  systems  [48],  but 
neither  is  currently  capable  of  supporting  heavy  or  awkwardly  shaped  targets.  The 
basic  measurement  setup  is  shown  in  Figure  1.1.  The  top  picture  is  an  aerial  view  of 
an  outdoor  range  in  New  Mexico.  A  bank  of  antennas  is  located  at  one  end  of  the 
range  (lower  right  in  Figure  1.1).  The  different  antennas  allow  for  different  frequency 
bands  to  be  measured.  The  target  is  located  opposite  the  antenna  bank  (lower  left 
in  Figure  1.1).  Note  the  absence  of  any  significant  structures  surrounding  the  pylon 
and  radar. 

The  target  can  be  rotated  and  inclined  to  allow  measurements  of  all  desired 
azimuth  and  elevation  angles.  The  pylon  itself  does  not  rotate,  but,  like  all  objects, 
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Figure  1.1:  RCS  Measurement  Setup  [4-6] 


it  scatters  incident  energy.  Pylons  have  been  shaped  to  enable  them  to  support 
significant  weight  while  minimizing  RCS  in  the  backscattering  direction.  Additionally, 
RAM  has  also  been  incorporated  into  the  pylons’  designs  to  help  reduce  the  scattered 
energy. 

The  collected  data  is  processed  to  calculate  the  DUT’s  RCS.  The  calculation 
used  to  determine  the  RCS  is  a  vector  background  subtraction  defined  as  [24] 


cr  = 


Ec  ~  Ecb 


aCAL, 


(1.1) 


s 

where  a  is  the  measured  RCS,  ET  is  the  scattered  field  when  the  target  is  mounted 

S  S 

on  the  pylon,  ETB  is  the  scattered  field  when  the  target  is  not  on  the  pylon,  Ec  is 
the  scattered  field  from  a  calibration  target,  ECB  is  the  scattered  field  when  the  cal¬ 
ibration  target  is  removed,  and  a  cal  is  the  calculated  RCS  of  the  calibration  target. 
Additionally,  the  calibration  target  is  typically  a  simple  shape  with  an  easily  theo¬ 
retically  determined  RCS.  This  calibration  is  done  to  identify  and  remove  sources  of 
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scattered  energy  other  than  the  DUT  as  well  as  compensate  for  systematic  errors  due 
to  the  non- ideal  radar.  However,  the  calibration,  while  good,  is  limited.  Equation  1.1 
cannot  be  used  to  calibrate  scattered  energy  resulting  from  a  target-pylon  interaction. 
This  is  when  the  scattered  energy  from  the  target  (pylon)  strikes  the  pylon  (target) 
and  results  in  a  signal  measured  by  the  radar.  Such  a  return  is  only  present  when 
the  target  is  mounted  on  the  pylon.  Thus,  it  cannot  be  corrected  via  a  calibration. 
Typically,  fields  from  this  type  of  interaction  are  small  compared  to  the  target’s  scat¬ 
tered  field.  However,  for  low  RCS  targets,  the  interaction  can  be  on  the  same  order  of 
magnitude.  Therefore,  a  way  to  reduce  the  bistatic  RCS  of  the  pylon  is  required.  A 
bistatic  reduction  is  necessary  because  the  energy  scattered  from  the  target  can  strike 
the  pylon  from  a  large  number  of  angles.  RAM  and  shaping  have  been  successfully 
used,  but  as  the  DUT  RCS  continues  to  decrease,  an  alternate  way  to  control  the 
component  of  the  scattered  field  resulting  from  a  target-pylon  interaction  is  required. 

1.2  Transformation  Optics 

Transformation  optics  is  a  relatively  new  field  that  provides  the  fundamental 
theory  enabling  precise  control  of  electromagnetic  waves.  Control  of  electromagnetic 
waves  is  certainly  not  a  new  technology.  Waveguides  and  fiber  optics  have  been  doing 
just  that  for  over  a  century.  A  key  distinction  is  waveguides  and  fiber  optics  are 
guiding  structures  operating  such  that  their  boundaries  confine  the  fields  within  a 
desired  space.  Transformation  optics  uses  a  smooth  variation  in  the  media  consti¬ 
tutive  parameters  to  steer  the  fields  in  a  desired  manner.  The  precision  with  which 
transformation  optics  allows  one  to  control  an  electromagnetic  field  is  unprecedented 
and  could  lend  itself  to  the  target-pylon  scattering  reduction  problem. 

Transformation  optics  works  because  geometric  rays  propagate  along  a  given 
trajectory  and  obey  Fermat’s  principle.  Fermat’s  principle  states  light  waves  of  a  given 
frequency  propagate  along  the  path  between  two  points  which  takes  the  least  time  [64], 
For  an  isotropic,  homogeneous  medium,  the  result  is  that  light  rays  propagate  in  a 
straight  line.  However,  when  the  medium  is  anisotropic  and  inhomogeneous,  the  path 
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which  the  rays  travel  can  be  quite  complex.  Thus,  by  controlling  the  constitutive 
material  parameters  (/Tr,  er ),  electromagnetic  energy  can  be  guided  in  any  way  one 
sees  fit.  But  how  does  one  exactly  design  a  desired  electromagnetic  response  in  a 
system?  Specifically,  other  than  using  trial  and  error,  how  can  one  design  the  required 
constitutive  parameters  of  a  medium  to  result  in  a  desired  electromagnetic  effect? 


Theoretically,  it  is  actually  quite  simple.  First,  one  develops  a  transformed  space 
in  which  electromagnetic  waves  propagate  in  a  desired  manner.  A  generic  transformed 
space  is  shown  in  Figure  1.2.  This  new  space  is  then  related  to  Cartesian  coordinates 


Figure  1.2:  Generic  transformation  space  [93] 


using  a  coordinate  transformation.  Ward  and  Pendry  showed  Maxwell’s  equations  are 
invariant  under  any  type  of  coordinate  transform  i.e.  the  equations  are  the  same  in  all 
coordinate  systems  with  only  the  permittivity  and  permeability  changing  values  [97]. 
Thus,  in  the  transformed  space,  Maxwell’s  equations  correctly  describe  the  behavior 
of  the  electromagnetic  waves.  One  may  use  the  invariance  of  Maxwell’s  equations  to 
derive  a  material  with  constitutive  parameters  defined  using  a  permittivity  and  per¬ 
meability  tensor.  When  this  material  is  placed  in  Cartesian  coordinates,  the  resulting 
held  behavior  accurately  mimics  the  held  behavior  in  the  transformed  space  (reference 
Appendix  A).  The  material  defined  by  the  permittivity  and  permeability  tensors  is 
what  creates  the  desired  electromagnetic  effects.  The  permittivity  and  permeability 
tensors  are  easily  calculable  once  the  coordinate  transformation  has  been  defined. 

The  transformed  space  can  encompass  any  type  of  electromagnetic  behavior  one 
desires  provided  a  one-to-one  transformation  exists  between  the  coordinate  systems. 
Simple  waveguide  bends,  held  concentrators,  and  space  that  contains  holes  where  no 
radiation  is  present  are  just  some  examples  that  have  been  simulated  in  computational 
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software  packages.  These  transformed  spaces  and  simulation  results  are  shown  in 
Figure  1.3. 


Figure  1.3:  Transformation  spaces  such  as  a  waveguide  bend  [93],  a  held  concen¬ 
trator  [73],  and  a  cylindrical  cloak  [56] 


At  this  point,  it  is  instructive  to  discuss  precisely  how  the  medium’s  constitu¬ 
tive  parameters  are  derived.  As  an  example,  consider  electromagnetic  cloaking.  The 
ability  to  cloak  an  object  using  metamaterials  (defined  in  Chapter  III)  was  first  dis¬ 
cussed  by  Pendry  et  al.  [72]  and  Leonhard t  [55]  in  2006.  The  techniques  discussed 
by  each  are  similar,  but  this  research  focuses  on  Pendry’s  method,  which  uses  the 
transformation  optics  approach  described  above.  An  electromagnetic  cloak  guides 
energy  around  a  particular  region  much  like  flowing  water  is  guided  around  a  stone. 
The  hidden  region  is  void  of  electromagnetic  energy,  meaning  an  object  can  be  placed 
in  the  hidden  region  without  perturbing  the  held.  What  follows  below  is  an  example 
of  the  transformation  optics  approach  used  to  derive  the  material  parameters  for  an 
infinitely  long  cylindrical  electromagnetic  cloak. 

1.2.1  Transformation  Optics  Cloaking  Example.  Per  transformation  optics, 
the  behavior  of  electromagnetic  waves  in  a  transformed  coordinate  system  can  be 
modeled  in  Cartesian  coordinates  using  a  material  with  specihc  permittivity  and 
permeability  tensors  [72,97].  This  derivation  is  shown  in  detail  in  Appendix  A,  and 
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the  resulting  constitutive  parameter  tensors  are  given  by 


el]  =  gtJ  \ui  ■  («2  x  u3)\QlQ2Q3(QiQj)  \  (1.2) 


filJ  =  gv  \iii  ■  (u2  x  u3)\ QiQ2Qz{QiQj)  L,  (1.3) 


where  glJ  are  the  components  of  the  inverse  of  the  coordinate  system’s 
which  is  dehned  as 


f  Ui  ■  Ui  U\  ■  U2  Ui  ■  u3  ^ 


metric  tensor 


9  = 


u2  ■  Ui  u2  ■  U2  u2  ■  u3 


(1.4) 


\  U3  ■  Ml  U3  U2  u3-u3  J 

Note  Ui  are  the  unit  vectors  in  the  i  —  1,  2,  3  direction  in  the  transformed  coordinate 
system,  and 


Qij 


dx  dx  dy  dy  dz  dz 
dqi  dqj  +  dq{  dq3  +  dqt  dq3  ’ 


Qi  Qiii  Qi  \J  Qii* 


(1.5) 

(1.6) 


An  electromagnetic  cloak  can  be  developed  by  creating  a  transformed  coordinate 
system  that  contains  voids  where  electromagnetic  energy  will  not  propagate.  The 
electromagnetic  field  behavior  in  the  transformed  coordinate  system  is  then  mimicked 
in  Cartesian  coordinates  using  a  material  with  permeability  and  permittivity  whose 
properties  are  described  by  Equations  1.2  and  1.3. 

Consider  a  transformed  cylindrical  coordinate  system  with  coordinates  ( r 6',  z') 
such  that  all  points  in  Cartesian  space  where  r  <  b  are  mapped  to  the  annular  region, 
a  <  r'  <  b.  This  can  be  written  mathematically  as 


r  =  ^1  -  r  +  a,  (1.7) 

where  r'  is  the  radial  location  in  the  transformed  coordinate  system  and  r  is  the  radial 
location  in  a  Cartesian  coordinate  system.  The  result  is  a  transformed  coordinate 
system  where  there  are  no  points  in  the  region  r'  <  a.  This  is  shown  in  Figure  1.4.  In 


r  =  b 


r'=  b 


Cartesian  Coodinates 


<^> 


Transformed  Coodinates 


Figure  1.4:  Transformed  Coordinate  System 


the  transformed  coordinate  system,  no  electromagnetic  energy  will  propagate  in  the 
region  r'  <  a  because  this  region  theoretically  does  not  exist.  Space  is  curved  around 
it. 

The  permittivity  and  permeability  tensors  which  electromagnetically  mimic  the 
curvature  of  the  transformed  space  can  be  found  as  follows.  The  mapping  from  the 
transformed  coordinate  system  where  r'  <  a  does  not  exist  to  Cartesian  coordinates 
can  be  written  as 


x  = 


a)b 


cos  O' , 


V  = 


a)b 


sin  O' , 


z  =  z 


(1.8) 


Note  the  transformed  coordinate  system  is  an  orthogonal  coordinate  system  with  unit 
vectors  f ,  0,  and  z.  Thus,  the  metric  tensor  is 


r  ■  r  r  ■  6  r  ■  z 

*->■ 

o 

o 

9  = 

0 ■ f  9-9  9 ■ z 

= 

0  1  0 

,  z  ■  r  z  ■  9  z  ■  z  , 

o  0  1 

(1.9) 


This  result  simplifies  the  expressions  in  Equations  1.2  and  1.3  due  to  the  fact 


f  ■  (  6  x  z 


=  5, 


ij, 


(1.10) 
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where  5,-j  is  the  Kronecker  delta  function.  Next,  the  Q\  values  can  be  found. 


Q 


2 

1 


/  rv  \  2  /  r\  \  2  /  rv  \  2 

/  ox\  /  oy  \  +(oz\ 
\8r  J  \8r  J  \8r  J 


(l.n) 


o^U-Y+U^V+U-Y 

\  r  89  )  \  r  89  J  \  r  89  J 

Q3  =  1 


The  partial  derivatives  can  be  expressed  as 


8x 

dr 


b 

b  —  a 


cos  6, 


8y_ 

dr 


b 

b  —  a 


sind, 


8z 

dr 


(1.12) 

(1.13) 


(1.14) 


1 8x 
r  89 


(r  —  a)b  sin  9  1  8y  (r  —  a)b  cos  9  1  8z 

b  —  a  r  ’  r  89  b  —  a  r  ’  r  89 


8x  8y  n  dx 

~8~z  =  ’  ~8z  =  ’  ~8~z 


Multiplying  out,  the  result  is 


Q  i  — 


b 

b  —  a 


Q2 


r  —  a  b 
r  b  —  a 

Q3  —  1- 


(1.15) 

(1.16) 

(1.17) 

(1.18) 
(1.19) 


Using  these  values  in  Equations  1.2  and  1.3  results  in  the  following  for  the  permittivity 
and  permeability  tensors  [79]. 


r  —  a 

£  y  (Jjy 

r 


(1.20) 


r 

£e  =  A4#  = - 

r  —  a 


r  —  a 
r 


b 

b - 


(1.21) 

(1.22) 
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These  can  be  rewritten  as  a  tensor  where  the  non-diagonal  terms  in  the  tensor  matrix 
are  zero. 

(  r~a  Q 
a 

jl  =  £=  0 

r— i 

v  o  o 

A  material  with  the  parameters  shown  in  Equation  1.23  and  immersed  in  free  space 
will  guide  all  electromagnetic  energy  around  the  region,  r  <  a,  much  like  the  energy 
would  propagate  in  the  transformed  coordinate  system  shown  in  Figure  1.4.  Thus, 
any  object  placed  in  this  region  will  not  scatter  any  electromagnetic  energy.  Hence, 
the  material  parameters  shown  in  Equation  1.23  effectively  define  an  electromagnetic 
cloak. 


(  6  j  ^  r— a  I 
\ b—a)  r  ) 


(1.23) 


1.2.2  Cloaking  and  Transformation  Optics.  As  shown  in  the  previous  sec¬ 
tion,  it  is  theoretically  possible  to  guide  incident  electromagnetic  energy  around  an 
object  such  that  the  object  has  no  scattered  field.  One  simply  needs  to  use  a  material 
having  the  constitutive  properties  described  by  Equation  1.23.  Transformation  optics 
led  to  this  result  and,  in  essence,  showed  cloaking  is  theoretically  possible  [72],  How¬ 
ever,  this  does  not  mean  the  hard  part  is  done.  Quite  the  contrary,  the  difficulty  lies  in 
developing  a  material  with  the  desired  constitutive  parameters.  Note  the  spatial  vari¬ 
ation  in  the  cylindrical  cloak’s  material  parameters.  Adding  to  the  complexity  is  the 
material  parameter  anisotropy.  A  material  with  the  properties  shown  in  Equation  1.23 
does  not  exist  naturally.  Fortunately,  advances  in  micro-  and  nano-fabrication  meth¬ 
ods  have  allowed  the  creation  of  man-made  materials  using  sub- wavelength  structures 
with  the  desired  material  properties  dictated  by  transformation  optics  [83].  Such  ma¬ 
terials  are  commonly  referred  to  as  metamaterials.  Metamaterials  are  the  enabling 
building  blocks  to  a  number  of  applications  spawned  from  transformation  optics. 
Metamaterials  will  be  discussed  in  detail  in  Chapter  III. 

Metamaterials  do  not  yet  enable  one  to  manufacture  an  ideal  cloak  with  the  pa¬ 
rameters  describe  by  Equation  1.23.  At  r  close  to  a,  the  diagonal  terms  in  Equation 
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1.23  are  going  to  either  zero  or  infinity.  A  material  with  infinite  permittivity /per¬ 
meability  will  likely  never  be  possible  to  manufacture.  Also,  all  existing  materials 
have  at  least  some  type  of  loss.  Creating  lossless  materials,  particularly  those  with 
magnetic  effects,  is  extremely  difficult.  Additionally,  the  required  material  param¬ 
eters  have  anisotropic  and  spatially  varying  /i  and  e,  which  is  quite  challenging  to 
make.  However,  simplifications  to  the  material  parameter  set  can  be  made  (Section 
2.2).  These  simplifications  result  in  less-than-ideal  cloaking  performance,  but  the 
end  result  does  maintain  some  of  the  ideal  cloak’s  electromagnetic  wave-controlling 
properties.  A  simplified  cylindrical  cloak  with  a  material  parameter  set  derived  from 
the  ideal  parameters  shown  in  Equation  1.23  has  recently  been  manufactured  and 
tested  [79]  with  promising  results. 

1.2.3  Cloaking  and  the  Speed  of  Light.  It  has  been  shown  it  is  possible  to 
cloak  a  region  of  space  such  that  an  observer  would  not  see  any  difference  in  the 
electromagnetic  fields  when  an  object  is  placed  in  this  hidden  region.  This  seems  to 
violate  the  fact  that  energy  cannot  propagate  faster  than  the  speed  of  light.  After 
all,  the  energy  must  be  bent  around  an  object  and  maintain  the  same  relative  phase 
as  the  energy  propagating  in  free  space.  Since  curving  around  an  object  requires  the 
energy  to  travel  a  further  distance,  it  seems  that  the  energy  must  propagate  faster 
than  the  speed  of  light.  However,  this  is  not  exactly  how  the  cloaking  process  works. 
The  energy  does  have  to  travel  a  longer  distance.  However,  the  cloak  is  not  trans¬ 
porting  energy  faster  than  the  speed  of  light.  Rather,  stored  energy  built  up  during 
the  transition  from  the  transient  to  the  steady-state  phase  allows  only  one  specific 
frequency’s  phase  fronts  to  exceed  the  speed  of  light  [2],  This  was  demonstrated  by 
Liang  et  al.  using  a  finite-difference  time-domain  (FDTD),  with  the  results  shown  in 
Figure  1.5.  Note  how  it  does  take  some  time  for  the  cloak  to  reach  its  stable  state. 
The  time  in  Figure  1.5  image  (a),  where  the  incident  wave  first  reaches  the  cloak  until 
steady-state  is  reached  in  image  (f)  is  approximately  15  periods  of  the  the  incident 
field  [59], 
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Figure  1.5:  Cylindrical  cloak  response  in  the  transient  to  steady-state  phase.  [59] 

While  still  a  long  way  from  implementing  a  cloaking  device  as  seen  in  Star  Trek, 
cloaks  present  a  new  paradigm  in  terms  of  RCS  reduction,  whether  it  be  creating  an 
ideal  cloak  (not  likely  in  the  near  future)  or  using  a  modified  cloaking  structure  in 
conjunction  with  shaping  and  RAM  to  further  reduce  an  object’s  overall  signature. 
Developing  a  cloaking  mechanism  for  a  support  pylon  on  an  RCS  range  could  help 
reduce  the  target-pylon  interactions  which  result  in  the  undesired  scattered  fields  due 
to  target-pylon  interactions  discussed  in  Section  1.1. 

1.3  Summary  of  Research  Goals 

The  intent  of  this  research  is  to  investigate  whether  cloaks  are  a  viable  option 
for  the  stated  long-term  problem.  Obviously,  a  three-dimensional  cloak  would  be  re¬ 
quired  for  any  real  implementation.  For  this  research,  however,  only  two-dimensional 
cylindrical  cloaks  are  considered.  This  is  due  to  the  fact  computer  requirements  for 
three-dimensional  cloak  simulations  are  rather  extensive,  whereas  two-dimensional 
simulations  are  easily  performed  on  a  standard  desktop  computer.  The  work  is  likely 
extendable  to  the  three-dimensional  case,  although  there  are  definite  issues  which 
must  be  considered  (Section  7.2).  Additionally,  a  unique  way  to  design  metamateri¬ 
als  to  increase  their  bandwidth  is  investigated. 
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There  are  two  main  thrusts  of  this  research.  First,  a  way  to  develop  simplified 
material  parameter  sets  for  cylindrical  cloaks  is  investigated.  This  is  important  be¬ 
cause  for  the  parameter  set  defined  by  Equation  1.23,  at  the  location  r  =  a,  /ir ,  er,  /iZl 
and  ez  all  equal  zero  while  [iq  and  £q  are  both  infinite.  These  values  are  unattainable 
no  matter  how  evolved  metamaterial  manufacturing  capability  becomes.  There  have 
been  some  generic  simplified  parameter  sets  published  in  the  literature.  These  simpli¬ 
fied  parameter  sets  maintain  some  cloaking  capability,  but  because  their  constitutive 
parameters  are  not  ideal,  cloaking  functionality  is  degraded  (Section  2.2).  A  method 
to  define  a  simplified  parameter  set  based  on  the  existing  metamaterial  manufacturing 
capabilities  has  not  been  developed.  Such  a  process  would  enable  cylindrical  cloak 
parameters  to  be  defined  in  terms  of  what  is  achievable,  thereby  not  putting  limits  or 
requirements  on  the  manufacturing  processes.  As  the  ability  to  manufacture  meta- 
materials  continues  to  advance,  material  parameter  sets  with  more  difficult  values 
will  be  able  to  be  obtained. 

The  second  thrust  of  this  research  involves  increasing  the  effective  bandwidth 
of  cloaks.  Note  the  ideal  cylindrical  cloak’s  material  parameters  shown  in  Equation 
1.23  are  independent  of  frequency.  Hence,  in  theory  the  cloak  is  wide  band  and 
would  be  well  suited  for  helping  to  reduce  a  pylon-target  interaction.  However,  cur¬ 
rent  research  has  shown  passive  metamaterials  used  to  realize  a  cloak  have  a  very 
narrow  operational  bandwidth  (Chapter  III).  Therefore,  a  cloak  constructed  using 
these  metamaterials  would  be  operationally  limited  to  a  small  range  of  frequencies. 
The  narrowband  nature  of  metamaterials  does  create  a  problem  because  RCS  ranges 
operate  over  a  significant  bandwidth.  A  narrow  band  solution  would  not  be  of  much 
use.  This  research  investigates  a  unique  way  to  increase  the  bandwidth  of  a  meta- 
material.  Making  the  building  blocks  have  a  broadband  response  would  result  in  the 
cloak  being  operational  over  a  larger  band  of  frequencies  and  would  help  make  cloaks 
a  more  viable  option  for  reducing  RCS  measurement  error. 
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1.4  Dissertation  Organization 

This  dissertation  is  organized  as  follows.  Background  information  on  cloaking 
theory,  the  cylindrical  cloak,  simplified  cloaks,  and  alternatives  to  cloaking  are  covered 
in  Chapter  II.  Chapter  III  examines  fundamentals  in  metamaterials  and  how  they 
are  designed  in  order  to  create  artificial  magnetic  and  electric  effects.  Additionally, 
common  methods  used  to  measure  the  constitutive  properties  of  metamaterials  are 
explained.  Chapter  IV  derives  a  new  constraint  equation  on  the  material  parameters 
for  ideal  cylindrical  cloaks,  which  is  then  used  as  the  foundation  to  develop  simpli¬ 
fied  material  parameter  sets.  Chapter  V  shows  how  a  Green’s  function  formulation 
can  be  used  to  decrease  solution  time  for  a  cylindrically  cloaked  perfect  electrically 
conducting  (PEC)  cylinder.  Chapter  VI  investigates  a  method  to  design  and  possibly 
increase  the  bandwidth  of  metamaterials.  Finally,  in  Chapter  VII,  conclusions  for  the 
research  are  summarized  with  recommendations  for  future  research. 
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II.  Cloaking  Background 

The  transformation  optics  design  approach  discussed  in  Chapter  I  provides  a  recipe 
for  various  types  of  electromagnetic  field  control.  This  research  focuses  on  electro¬ 
magnetic  cloaking,  which  was  first  put  forth  by  Pendry  et  al.  in  2006  [72],  Since  the 
publication  of  this  ground-breaking  work,  there  have  been  numerous  papers  published 
analyzing  the  behavior  of  ideal  cloaks  using  common  electromagnetic  analysis  tech¬ 
niques.  The  relevant  works  are  discussed  below.  As  noted  in  Section  1.2.2,  materials 
for  an  ideal  cylindrical  cloak  do  not  exist,  neither  naturally  nor  can  they  be  perfectly 
manufactured.  This  limitation  necessitates  simplified  parameter  sets.  A  number  of 
simplified  parameter  sets  for  a  two-dimensional  cylindrical  cloak  have  been  devel¬ 
oped.  These  are  discussed  in  addition  to  their  various  short-comings.  Finally,  other 
cloaking  options  not  based  on  transformation  optics  are  briefly  discussed  and  docu¬ 
mented.  The  limitations  associated  with  ideal  cloaking  that  involve  the  design  and 
manufacture  of  metamaterials  will  be  discussed  in  Chapter  III. 

2.1  Perfect  Cloaking  Theoretical  Analysis 

Since  Pendry  et  al.’s  initial  paper  in  2006,  there  has  been  a  significant  effort 
confirming  that  perfect  cloaking  is  theoretically  possible,  assuming  the  ideal  consti¬ 
tutive  parameters  dictated  by  transformation  optics  could  be  achieved.  Schurig  et 
al.  developed  a  method  to  perform  ray-tracing  within  a  cloak  in  order  to  confirm  the 
cloak  behaves  as  theoretically  derived  [81].  For  spherical  and  cylindrical  cloaks,  they 
showed  via  ray-tracing  the  complex  material  acts  as  a  perfect  cloaking  mechanism 
for  the  desired  hidden  region  while  resulting  in  no  perturbation  to  the  incident  ray 
trajectory  outside  the  cloaking  body. 

Leonhardt  and  Philbin  demonstrated  how  transformation  optics  and  the  associ¬ 
ated  behavior  of  the  electromagnetic  fields  can  be  described  using  the  general  theory  of 
relativity  [56].  They  developed  a  formulation  which  takes  a  desired  function,  whether 
it  be  cloaking,  perfect  lenses,  or  the  behavior  of  artificial  black  holes,  and  finds  the 
properties  of  the  material  needed  to  generate  the  desired  behavior.  They  showed  the 
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behavior  of  electromagnetic  fields  in  cloaks  can  be  included  and  described  within  the 
framework  of  general  relativity,  which  further  solidified  the  validity  of  electromagnetic 
cloaking  as  described  by  transformation  optics. 

Chen  et  al.  performed  a  full  wave  Mie  scattering  analysis  on  a  spherical  cloak 
[18].  They  quantitatively  solved  for  the  scattered  fields  from  an  ideal  spherical  cloak 
and  determined  the  scattering  would  be  identically  zero.  When  loss  was  introduced  to 
the  cloaking  material,  bistatic  scattering  resulted,  with  larger  losses  equating  to  larger 
scattered  fields.  However,  the  addition  of  loss  did  not  affect  the  back-scattered  field 
in  that  the  monostatic  return  was  still  zero.  This  result  is  very  different  from  that 
of  regular  particles  and  applies  only  to  the  spherical  cloak.  When  introducing  loss  to 
the  ideal  parameters  for  a  cylindrical  cloak,  the  monostatic  return  is  not  identically 
zero. 

Ruan  et  al.  used  cylindrical  wave  expansion  to  solve  for  the  scattered  field  from 
an  ideal  cylindrical  cloak.  They  also  solved  for  the  held  transmitted  into  the  hidden 
region.  They  confirmed  by  applying  boundary  conditions  the  ideal  cloak  is  perfect  by 
proving  the  coefficients  for  the  scattered  held  from  the  cloak  and  the  transmitted  held 
into  the  hidden  region  were  all  zero  [77] .  This  proved  a  cylindrical  cloak  with  the  ideal 
parameters  shown  in  Equation  1.23  has  no  rehected  held  in  addition  to  providing  a 
hidden  region  (r  <  a  in  Figure  1.4)  which  is  completely  shielded  from  electromagnetic 
energy. 

Weder  studied  hrst-order  and  higher-order  spherical  cloaks.  He  proved  for  any 
frequency  that  ideal  cloaks  have  no  scattered  held.  Additionally,  he  showed  that  no 
incident  energy  can  penetrate  into  a  cloak’s  hidden  region,  and  that  if  a  source  were 
placed  in  the  hidden  region,  its  energy  would  not  leave  the  concealed  area  [98].  This 
makes  sense  because  reciprocity  holds  for  spherical  cloaks  since  the  permittivity  and 
permeability  tensors  are  symmetric  [52], 

Zhang  et  al.  developed  the  equations  to  formulate  the  material  parameters 
necessary  to  cloak  an  object  in  a  slowly  varying,  multilayered,  inhomogeneous  envi- 
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ronment  [106].  The  work  by  Pendry  et  al.  assumed  a  cloak  immersed  in  homogeneous 
free  space.  Zhang  et  al.' s  analysis  is  similar  to  the  original  work  by  Ward  and  Pendry. 
Their  simulations  showed  successful  cloaking  of  an  object.  The  application  for  their 
work  is  cloaking  objects  in  a  layered  media  or  at  an  interface  between  two  media, 
such  as  a  ship  at  sea  or  a  building  in  the  hot  desert  air. 

A  computer  simulation  of  the  perfect  two-dimensional  cylindrical  cloak  was  done 
by  Cummer  et  al.  They  performed  a  full-wave  finite  element  method  (FEM)  simu¬ 
lation  of  the  ideal  two-dimensional  cylindrical  cloak  using  the  Cornsol  Multiphysics 
FEM-based  electromagnetics  solver  [28].  They  simulated  the  two-dimensional  cylin¬ 
drical  cloak  for  both  lossless  and  lossy  materials  using  a  transverse  magnetic  (TMZ) 
incident  wave  (an  electromagnetic  wave  with  only  a  ^-component  for  the  electric  field 
vector).  Their  simulation  results  clearly  showed  the  cloak  operating  as  theoretically 
predicted  with  some  degradation  in  performance  when  loss  was  introduced. 

There  are  a  number  of  papers  which  derive  the  theoretical  equations  for  various 
cloaking  geometries.  Ma  et  al.  [60]  used  the  transformation  optics  algorithm  described 
by  Pendry  et  al.  to  derive  the  material  parameter  equations  for  an  elliptical  cylindrical 
cloak  with  similar  results  shown  in  [45].  Kwon  and  Werner  did  a  similar  analysis, 
but  they  considered  an  eccentric  elliptic  electromagnetic  cloak  [53].  Rahm  et  al. 
designed  and  simulated  a  square  cloak  and  a  cylindrical  concentrator,  which,  instead 
of  cloaking  a  certain  region,  focuses  fields  from  one  region  into  another  [73].  Jiang  et 
al.  considered  conformal,  arbitrarily  shaped  cloaks  [44],  These  papers  all  performed 
simulations  using  the  Cornsol  Multiphysics  software  package,  and  the  results  clearly 
showed  the  cloaks  (or  the  concentrator)  working  as  predicted  by  the  original  theory. 

Zhao  et  al.  performed  a  full-wave  FDTD  analysis  of  a  two-dimensional  cylin¬ 
drical  cloak.  They  used  the  Drude  dispersion  model  to  represent  the  permittivity 
and  permeability  of  the  cloak’s  material  parameters.  As  with  other  simulations,  they 
found  a  cloak  with  ideal  parameters  effectively  hides  an  object  placed  within  the 
cloaking  shell  from  incident  electromagnetic  energy.  Liang  et  al.  also  performed 
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FDTD  simulations  on  a  cylindrical  cloak  and  found  the  cloak  to  work  as  expected. 
However  they  noted  a  strong  forward  scattering  from  the  cloak  during  the  dynamic 
processes  when  the  incident  waveform  first  strikes  the  cloak,  an  effect  that  can  be 
controlled  by  varying  the  dispersive  parameters  in  the  Drude  model  [59]. 

These  theoretical  results  support  the  original  derivation  that  a  cloak  of  a  specific 
geometry  with  material  parameters  derived  according  to  the  transformation  optics 
method  put  forth  by  Pendry  et  al.  does  indeed  result  in  a  perfectly  cloaked  region. 

2.2  Simplified  Cylindrical  Cloaks 

This  section  focuses  solely  on  two-dimensional  cylindrical  cloaks.  The  units  for 
this  geometry  in  this  work  are  ( r,(j),z ),  which  is  consistent  with  the  work  published 
by  Schurig  et  al.  [79].  The  cylindrical  cloak  has  a  hidden  region  located  at  r  <  a, 
where  a  is  the  inner  boundary.  Objects  placed  in  the  hidden  region  are  completely 
shielded  from  electromagnetic  energy.  The  outer  boundary  of  the  cloak  is  located  at 
r  =  b.  Additionally,  all  analysis  in  this  work  assumes  plane  wave  incidence.  This 
geometry  is  shown  in  Figure  2.1. 


Plane  Wave 
Incidence 

0i=  180° 


Figure  2.1:  Two-dimensional  cylindrical  cloak  geometry 
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When  examining  the  material  properties  of  a  cylindrical  cloak  (Equation  1.23), 
it  is  obvious  such  a  material  would  not  likely  occur  naturally.  First  note  the  cloak 
is  lossless.  All  materials  have  loss,  and  it  is  very  difficult  to  manufacture  magnetic 
materials  with  even  small  losses.  Also,  (j,r,  er,  Hz  and  all  equal  zero  at  r  =  a. 
Additionally  He  and  Eg  are  each  infinite  at  r  =  a.  Thus,  not  only  are  these  materials 
not  naturally  occurring,  it  is  not  currently  possible  to  manufacture  them,  nor  is  it  very 
likely  technology  would  ever  enable  the  manufacture  of  infinite  values.  Because  of  this, 
simplified  parameter  sets  have  been  derived  with  the  intent  of  creating  manufacturable 
constitutive  parameter  values  while  limiting  the  reduction  in  cloak  functionality  as 
much  as  possible. 

As  a  way  to  reduce  the  number  of  constitutive  parameters  required  to  realize 
a  manufacturable  cloak,  the  incident  field  can  be  decomposed  into  transverse  electric 
(TE)  and  transverse  magnetic  (TM)  field  components.  Thus,  for  TE2  fields  only  Hz, 
er,  and  Eg  are  required  when  analyzing  field  behavior.  TM2  fields  require  ez,  Hr ,  and 
He- 

The  first  sets  of  simplified  material  parameters  for  a  cylindrical  cloak  were 
developed  for  specific  incident  field  polarizations  with  the  goal  of  satisfying  the  same 
governing  wave  equation  within  the  simplified  cloak  that  is  satisfied  in  the  ideal 
cloak  [79].  For  an  assumed  incident  field  type,  Maxwell’s  equations  can  be  used  to 
define  a  wave  equation  that  governs  the  field  behavior  within  a  given  space.  Assuming 
TM2  incidence,  Maxwell’s  equations  can  be  expressed  as 


Ek  = 


JL0£z£or 


Hr  = 


d  (tHq 

dr 


dHr 

~w 


dEz 


Ho  = 


juHrHo'r  96 

1  dEz 


(2.1) 

(2.2) 

(2.3) 


juJHeHo  dr 

where  Ez  is  the  ^-component  of  the  electric  field,  Hr  is  the  r-component  of  the 
magnetic  field,  Hg  is  the  ^-component  of  the  magnetic  held,  r  is  the  radial  location, 
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e0  and  /iG  are  the  permittivity  and  permeability  of  free  space,  £z,  fie-,  and  jir  are  the 
relative  permittivity  and  permeability  tensor  values,  and  tv  is  the  angular  frequency  of 
the  incident  electromagnetic  held.  The  general  wave  equation  governing  the  behavior 
of  TM~  fields  within  a  complex  anisotropic  material  with  spatially  varying  constitutive 
parameters  can  be  developed  by  substituting  Equations  2.2  and  2.3  into  Equation  2.1. 


Ez  jueze0r  dr  \junoii0  dr 


jujUrHo'r  d6 


Note  that  j,  tv,  Ho,  and  £0  are  all  independent  of  r  and  6.  Therefore,  they  can  be 
removed  from  the  differentiation  operations. 


tv2ii0£0£zr  dr  \jio  dr 


r  dE, 


d  f  1  dEz 
86  \  /irr  dO 


Since  r  is  not  a  function  of  0,  it  can  be  removed  from  the  operation.  Additionally, 
k0  is  the  free-space  wave  number  defined  as  kQ  =  tv^/j,0£0.  Therefore,  Equation  2.5 
can  be  rewritten  as 


1  r  d  (  r  dEz\  1  1  d  (1  dEz\  2 

£zr  dr  \/ig  dr  )  +  £zr2  dO  \/ir  dO  )  ° 


(2.6) 


Equation  2.6  is  the  general  wave  equation  that  describes  the  behavior  of  a  TM2 
electromagnetic  held  in  a  cylindrical  anisotropic  media.  This  equation  will  be  used 
in  the  development  of  simplified  constitutive  parameter  sets  for  cylindrical  cloaks. 

Schurig  et  al.  were  the  hrst  to  derive  a  set  of  simplified  material  parameters  for 
a  two-dimensional  cylindrical  cloak  [79].  However,  when  deriving  the  wave  equation 
that  governs  the  helds  within  the  cloak,  the  procedure  they  used  assumed  a  priori 
He  was  constant.  Their  intent  was  to  eventually  simplify  fie  to  a  constant  value,  but 
doing  so  when  developing  the  wave  equation  was  mathematically  incorrect  [101,102], 
Because  of  this  error,  Schurig  et  al.  simply  removed  He  from  the  differentiation 
operation  with  respect  to  r  in  Equation  2.5.  Thus,  the  following  was  thought  to  be 
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the  wave  equation  for  TM2  fields  in  a  complex  anisotropic  medium. 


1  d2Ez  11  dEz  11  d2Ez  ,  9  „ 

- TT"  ^ - 7) - t - 2  ou2  “I”  —  0 

e2/xe  drz  £z^e  r  or  £zfxr  rz  dt)z 


(2.7) 


Note  Equation  2.7  does  not  equal  Equation  2.6  if  he  is  r-dependent.  Not  realizing 
this,  Schurig  et  al.  then  substituted  the  ideal  cylindrical  cloak’s  material  parameters 
for  a  TM2  incident  field  into  Equation  2.7.  The  following  is  the  result. 


b  —  a\2  d2Ez  fb-a\2ldEz  fb_-a\2  f  1  V  d2Ez 
b  )  dr 2  \  b  )  r  dr  \  b  )  \r  —  a)  d9'2  ° 


0  (2.8) 


Schurig  et  al.  believed  Equation  2.8  was  the  correct  wave  equation  for  TM2  fields  in  an 
ideal  cylindrical  cloak.  Their  goal  was  to  develop  a  simplified  cylindrical  cloak  whose 
internal  held  behavior  would  match  the  held  behavior  in  a  cloak  with  ideal  parameters. 
Therefore,  they  compared  Equations  2.7  and  2.8  and  concluded  the  following  were 
the  only  material  constraints  on  a  simplified  cylindrical  cloak’s  material  parameters 
for  TM2  incident  fields  in  order  for  the  electric  held  to  satisfy  the  same  wave  equation 
as  that  of  an  ideal  cloak  [79]. 


1 


(2.9) 


1  /  b  —  a  \  2  /  r 

£z.l-h  V  b  J  \r  -a 

By  examining  Equations  2.9  and  2.10,  Schurig  et  al.  developed  the  following  set  of 
material  parameters  for  a  simplified  cylindrical  cloak. 


(2.10) 


hr  — 


r  —  a 
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) 


he  —  1; 


(2.11) 


Note  that  he  is  indeed  a  constant  like  they  assumed,  but  because  this  was  assumed 
before  all  mathematical  operations  were  completed,  their  results  were  not  entirely 
correct.  The  material  parameters  shown  in  Equation  2.11  do  satisfy  Equations  2.9 
and  2.10,  and,  as  will  be  shown  in  Chapter  IV,  Equations  2.9  and  2.10  are  some  (but 
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not  all)  of  the  constraint  equations  on  the  material  parameters  for  an  ideal  cylindrical 
cloak.  Thus,  the  simplified  parameters  developed  by  Schurig  et  al. ,  though  obtained 
using  a  method  with  a  mathematical  error,  did  show  good  results.  Simulations  showed 
a  cloak  with  these  simplified  parameters  maintained  much  of  the  power-flow  bending 
and  low-reflection  characteristics  of  the  ideal  cloak  [28]. 

These  simplified  parameters  are  not  nearly  as  complex  as  the  ideal  material 
parameters  from  a  manufacturing  perspective,  as  is  now  position  invariant.  Addi¬ 
tionally,  Schurig  et  al.  defined  no  —  1  because  this  is  an  easily  obtainable  value.  The 
only  spatially  varying  parameter  is  /ir .  No  parameters  are  infinite  and  only  /ir  =  0  at 
r  =  a.  However,  in  order  to  realize  such  a  device,  one  further  simplification  had  to  be 
made.  Note  that  fir  is  still  spatially  varying.  Concentric,  homogeneous  layers  with 
Hr  =  nn  were  used,  where  fin  is  a  constant  value  for  the  nth  layer.  A  cloak  with  these 
simplified  parameters  was  manufactured  using  metamaterials  [79].  The  results  were 
not  as  good  as  the  simulated  values  with  varying  /ir ,  but  they  did  show  the  ability  to 
partially  cloak  an  object. 

The  simplified  parameter  set  shown  in  Equation  2.11  was  used  in  the  design  and 
simulation  of  a  cloak  manufactured  using  metamaterials  consisting  of  high  permittiv¬ 
ity  ferroelectric  rods.  Galliot  et  al.  used  metamaterial  building  blocks  consisting  of 
Ba^Sri-zTiOs  rods,  and  by  adjusting  the  rod  radii,  they  could  control  the  resonant 
frequency.  Much  like  Schurig  et  al.,  Galliot  et  al.  created  the  radial  variation  in  /ir 
by  creating  layers  of  concentric  rings.  The  difference  is  the  operating  frequency  of 
their  cloak  was  0.58  THz  compared  to  8.5  GHz  for  the  Schurig  et  al.  cloak.  Galliot 
et  al.  also  performed  simulations  using  the  commercial  FEM  software  package,  High 
Frequency  Structure  Simulator.  Their  work  was  unique  because,  unlike  Cummer  et 
al.  in  [28],  Gaillot  et  al.’s  simulated  the  individual  building  blocks  of  their  cloak  i.e. 
the  metamaterial  structures.  Cummer  et  al.  used  continuous  subdomains  in  their 
simulations.  The  results  of  Gaillot  et  aV s  work  showed  a  simplified  cloak  with  some 
of  the  ideal  characteristics  operating  in  the  THz  region  [35]. 
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The  original  simplified  parameter  set  was  used  to  develop  a  material  parameter 
set  for  construction  of  a  non-magnetic  simplified  cylindrical  cloak  for  TEZ  incident 
fields  at  optical  frequencies  [14].  A  non-magnetic  cloak  was  desired  due  to  the  dif¬ 
ficulty  in  manufacturing  materials  with  a  magnetic  response  at  optical  wavelengths. 
Recall  for  TEZ  incident  fields,  the  constitutive  parameters  required  are  )iz,  eg,  and  er. 
Following  a  derivation  process  similar  to  that  done  by  Schurig  et  al.  and  described 
above,  the  constraint  equations  for  TEZ  incident  fields  were  thought  to  be 
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(2.12) 
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/iz£r 


(2.13) 


Using  Equations  2.12  and  2.13  and  the  desire  to  limit  the  material  parameters  to 
non-magnetic  effects,  the  following  simplified  material  parameters  were  developed. 


Ah  1) 


(2.14) 


Much  like  the  results  in  [28] ,  simulations  showed  this  simplified  cloak  maintains  some 
of  the  characteristics  of  the  ideal  cloak  [14], 

There  have  been  improvements  to  the  original  simplified  parameters.  An  ob¬ 
vious  shortcoming  of  the  original  simplified  parameters  is  the  large  reflection  at  the 
cloak’s  outer  boundary.  The  ideal  cloak  has  an  impedance  matched  to  free  space  at 
r  =  b.  The  original  simplified  cloak  does  not,  resulting  in  a  significant  reflection  at 
the  cloaking  body  and  free  space  interface.  To  fix  this  problem,  it  was  noted  the 
transformation,  denoted  as  g{r'),  mapping  the  space  r'  <  b  to  the  cylindrical  shell 
a  <  r  <  b,  can  have  multiple  forms  [15].  The  linear  transformation  has  the  form 


r  = 


r'  +  a, 


(2.15) 
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while  a  quadratic  transformation  can  have  the  form 


r  =  g(r') 


b)  +  1 


r  +  a. 


(2.16) 


As  before,  a  and  b  are  the  inner  and  outer  radii  of  the  cloak  while  p  is  a  parameter 
which  will  be  determined  shortly.  It  can  be  shown  using  the  procedure  developed 
in  [72]  the  material  parameters  for  the  ideal  cylindrical  cloak  can  be  represented  as 
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dr'  '  fir 

A  reduced  parameter  set  was  desired  for  TM2  incident  fields.  When  simplifying  the 
ideal  parameter  set  shown  in  Equations  2.17  and  2.18,  the  constraints  defined  by 
Equations  2.9  and  2.10  were  used  as  the  only  limits  on  the  material  parameters. 
Hence,  by  setting  £z  —  1,  the  following  define  a  set  simplified  parameters. 
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To  match  the  impedance  to  free  space  at  r  =  b,  the  following  constraint  was  applied. 
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(2.20) 


Using  the  condition  imposed  in  Equation  2.20,  the  variable  p  in  Equation  2.16  is 
found  to  be  p  =  A.  For  g(r')  equal  to  that  shown  in  Equation  2.16  with  p  —  A, 
the  parameter  set  is  called  a  quadratic  cloak  [15].  The  quadratic  cloak  satisfies  the 
material  constraints  shown  in  Equations  2.9  and  2.10  and  has  a  matched  impedance 
at  r  =  b.  A  limit  on  the  quadratic  cloak  is  the  value  |  <  0.5  to  ensure  a  monotonic 
transformation  [15];  however,  its  performance  in  terms  of  reducing  the  scattering 
width  of  a  PEC  cylinder  was  better  than  that  of  the  cloak  with  the  original  simplified 
parameters. 
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A  second  set  of  material  parameters  with  a  matched  impedance  at  r  —  b  that 
satisfies  equations  (2.9)  and  (2.10)  was  developed  [102],  These  are  shown  below. 


/ir 


(2.21) 


To  the  author’s  knowledge,  the  parameter  set  in  Equation  2.21  was  not  derived  from 
a  governing  equation.  Rather,  the  values  seem  to  have  been  found  by  simply  using 
Equations  2.9  and  2.10  and  determining  values  which  satisfy  these  constraints  while 
also  having  a  matched  impedance  at  r  =  b.  For  this  work,  the  cloak  with  parameters 
shown  in  Equation  2.21  will  be  called  the  improved  cloak.  The  performance  in  terms 
of  reducing  the  scattering  width  of  a  cloaked  PEC  cylinder  by  the  improved  cloak 
and  the  quadratic  cloak  are  similar  for  certain  values  of  a  and  b.  However,  it  has 
been  shown  the  improved  cloak  has  a  more  consistent  performance  as  the  ratio  of  | 
varies  [102],  Additionally,  the  improved  cloak  has  only  one  spatially  varying  material 
parameter  while  the  quadratic  cloak  has  two,  making  the  improved  cloaks  parameters 
easier  to  manufacture  than  those  of  the  quadratic  cloak. 


2. 3  Cloaking  Limitations 

As  stated  in  the  previous  section,  the  ideal  cylindrical  cloak  has  constitutive 
parameters  equal  to  zero  or  infinity,  values  which  are  not  possible  to  obtain  and  are 
the  motivation  to  develop  simplified  parameter  sets.  It  is  obvious  since  the  ideal 
parameter  values  won’t  be  obtained,  there  will  be  a  reduction  in  cloak  performance. 
Isic  et  al.  analyzed  cloak  performance  based  on  the  inability  to  precisely  manufacture 
the  ideal  cloak’s  constitutive  parameters  [43].  Note  for  the  ideal  cloak,  gr  =  0  at  the 
inner  boundary,  r  =  a.  This  also  means  go  — >  oo  at  the  outer  boundary,  r  =  b.  It  is 
not  possible  to  achieve  these  values.  For  their  analysis,  Isic  et  al.  let  the  parameters 
for  a  TM2  incident  wave  be 

r  _  b(a-r1) 

gr  = - b-=^~  (2.22) 
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The  parameter,  r i,  is  a  measure  of  the  imperfectness  of  a  cloak.  An  ideal  cloak  has 
ri  =  0.  At  this  value  Equations  2.22  -  2.24  simplify  to  the  ideal  parameters.  Defining 
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and  since  n rideoi(r  =  a)  =  0,  they  solve  for  rq  to  be 


r  i  = 


ab5fir 


b  —  a(  1  —  8  nr) 


(2.26) 


For  a  TMZ  plane  wave  incident  on  the  imperfect  cloak,  and  with  a  PEC  cylinder  placed 
in  the  hidden  region,  they  solved  boundary  conditions  and  found  the  coefficients  for 
the  scattered  field.  For  small  ri,  the  dominant  scattering  term  is  the  zeroth  order 
mode,  which  matches  the  results  of  Ruan  et  al.  Isic  et  al.  defined  q  as  the  ratio  of 
the  scattering  width  of  an  uncloaked  PEC  cylinder  to  that  of  a  cloaked  PEC  cylinder 
and  find  that 


Q  = 


2.29  In  2(k0r1) 
n\n 


(2.27) 


If  one-order  of  magnitude  scattering  width  reduction  is  desired,  q  =  10.  The  param¬ 
eters  ri  and  5/j,r  can  be  found  using  Equations  2.26  and  2.27.  For  q  =  10,  A0  =  0.25 
m,  a  =  X0,  and  b  =  2X0,  they  found  that  Sfir  =  0.01.  This  means  in  order  to  reduce 
the  scattering  width  by  10  dB  using  a  cylindrical  cloak,  the  constitutive  parameter, 
Hr  =  0.01  at  r  —  a  [43].  Such  manufacturing  control  is  possible,  but  anything  more  is 
approaching  the  current  technological  limit.  In  [42],  Isic  et  al.  provide  a  more  mathe¬ 
matically  rigorous  derivation  of  the  scattering  from  imperfect  cloaks.  They  conclude 
that  a  PEC  object  placed  in  the  hidden  region  of  an  imperfect  cloak  has  a  reduction  in 
its  echo  width  by  a  factor  of  ^  [42],  They  consider  one  order  of  magnitude  reduction 
an  optimistic  result. 
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Ruan  et  al.  performed  a  similar  analysis  and  examined  the  effects  of  a  slight 
perturbation  in  the  location  of  the  cloak’s  inner  boundary.  They  found  even  small 
perturbations  result  in  scattering  from  the  cloaking  body.  A  perturbation  takes  the 
form  of  the  inner  cloak  boundary  being  located  at  a  distance  r  =  a  +  5  instead 
of  the  ideal  location  r  =  a.  Even  for  very  small  <5,  such  as  5  =  10_5a,  there  is  a 
noticeable  scattered  held  [77].  The  size  of  the  scattered  held  is  dependent  on  the 
cloak  geometry.  As  an  example,  consider  a  cloak  with  boundaries  a  =  A  and  b  =  2a. 
A  ^-directed  incident  and  scattered  electric  held  can  be  represented  using  cylindrical 
wave  expansion  as  [13] 

OO 

E'^E^aMk  r)e>"*,  (2.28) 


Ei  =  E0YJC„H^(kr)^, 


(2.29) 


where  E0  is  a  constant  value.  The  scattering  coefficient  for  each  order  can  be  defined 
as 

<  =  -■  (2-30) 

an 

For  5  =  10_5a,  | cvg |  =  0.175.  For  orders  where  n  ^  0,  \a^\  <  10^9,  meaning  the 
zeroth  order  is  the  dominant  scattering  term.  Additionally,  it  was  shown  as  6  — >  0, 
convergence  of  |ckq|  is  slow  [77]. 

In  addition  to  having  constitutive  parameters  which  are  unattainable,  perfect 
cloaks  also  have  a  bandwidth  issue.  Yao  et  al.  investigated  whether  or  not  the 
material  making  up  an  electromagnetic  cloak  could  be  frequency  invariant  [104],  They 
concluded  that,  due  to  causality,  the  cloak  must  be  dispersive.  A  nondispersive  cloak 
results  in  group  velocities  greater  than  the  speed  of  light.  Additionally,  they  found 
there  must  be  a  strong  absorbtion  at  the  cloak’s  frequency  of  operation.  This  results 
in  a  significant  forward  shadow.  They  did  conclude  the  cloak  is  an  effective  instrument 
to  reduce  backscatter  but  only  for  a  narrow  bandwidth. 


Chen  et  al.  had  similar  conclusions  regarding  the  bandwidth  as  a  result  of 
the  limit  imposed  on  group  velocities  [20].  They  considered  a  more  generic  form  of 
the  cylindrical  transformation  such  that  they  are  mapping  the  region  rQ  <  r  <  b  to 
a  <  r  <  b.  In  the  example  in  Section  1.2.1,  ra  =  0.  Chen  et  al.  showed  perfect 
cloaking  is  possible  only  at  one  single  frequency.  If  more  than  one  frequency  is 
considered,  the  group  velocities  exceed  the  speed  of  light.  Therefore,  they  considered 
a  cloak  such  that  r0  is  now  a  function  of  frequency,  rQ( u).  Note  rQ(u)  can  be  zero 
at  one  particular  frequency  and  also  have  small,  but  non-zero  values  throughout  a 
nearby  bandwidth.  They  simulated  a  cloak  with  the  following  material  parameters. 
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They  show  the  bandwidth  limitation  for  this  material  parameter  set  to  be 
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They  simulated  their  designed  cloak  and  found  it  to  have  an  operational  bandwidth 
from  8.5  -  8.75  GHz,  with  the  cross  section  of  a  PEC  cylinder  being  30%  compared 
to  that  of  an  uncloaked  cylinder. 

As  a  way  to  increase  the  operational  bandwidth  of  a  cylindrical  cloak,  Wang 
et  al.  proposed  creating  a  cylindrical  cloak  out  of  active  metamaterials  [96].  They 
designed  a  simplified  cylindrical  cloak  with  material  parameters  shown  in  Equation 
2.14  and  repeated  below  for  convenience. 


Note  for  their  effort,  they  assumed  a  TE2  incident  waveform.  For  their  cloak  im¬ 
plementation,  they  proposed  to  use  active  metamaterials.  Passive  metamaterials  are 
dispersive  (to  be  discussed  in  Chapter  III)  i.e.  their  constitutive  parameters  vary  as 
a  function  of  frequency.  For  the  geometry  of  their  cloak  and  assumed  incident  field 


29 


type,  the  gz  and  £e  values  will  be  determined  by  the  substrate  used  in  the  metamate¬ 
rial  manufacture,  and  the  substrate  is  not  dispersive  in  the  desired  frequency  range. 
However,  er  is  determined  by  the  configuration  of  the  metallic  structures  within  the 
metamaterial  and  is  highly  dispersive.  In  order  to  allow  a  band  of  operating  fre¬ 
quencies,  Wang  et  al.  proposed  using  a  variable  capacitor  mounted  between  metal 
strips.  The  capacitance  between  the  metal  strips  helps  control  er.  A  DC  bias  can  be 
used  to  change  the  capacitance  on  the  variable  capacitor.  As  frequency  is  changed, 
the  capacitance  is  changed  in  order  to  keep  the  er  at  the  required  value  dictated  by 
Equation  2.33.  Wang  et  al.  did  not  build  any  devices,  but  their  numeric  simulations 
showed  their  concept  is  valid.  The  instantaneous  bandwidth  of  the  cloak  was  not 
increased,  but  it  was  able  to  operate  over  a  wider  bandwidth. 

The  final  work  discussed  in  this  section  is  not  really  a  short-coming  of  cloaks, 
but  rather  a  unique  way  in  which  the  cloaking  capability  can  be  turned  off.  Chen  et 
al.  used  the  transformation  optics  approach  to  develop  a  material  that  when  placed 
in  the  hidden  region  of  an  ideal  cylindrical  cloak  results  in  the  object  within  the 
hidden  region  being  seen  by  incident  radiation.  A  PEC  cylinder  was  placed  in  the 
hidden  region  of  an  ideal  cloak  and  was  coated  with  a  specific  anisotropic  material 
with  a  negative  refractive  index.  Chen  et  al.  termed  this  material  the  anti-cloak. 
The  anti-cloak  effectively  annihilates  the  functionality  of  the  ideal  cloak  and  shifts 
the  hidden  region  out  to  the  cloak’s  boundary  thereby  making  it  visible.  They  proved 
the  existence  of  scattered  fields  by  matching  boundary  conditions  [21]. 

2-4  Alternate  Cloaking  Methods 

This  section  is  by  no  means  a  complete  listing  of  the  additional  work  going  on 
with  regards  to  cloaking,  but  it  does  provide  a  general  summary  of  additional  work 
being  done  in  this  field. 

The  work  done  by  Huang  et  al.  is  cloaking  as  described  by  transformation  optics. 
However,  to  implement  the  cloak,  they  did  not  use  anisotropic  materials.  Rather, 
they  assumed  a  TEZ  incident  wave  and  simulated  a  two-dimensional  cylindrical  cloak 
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realized  using  homogeneous  isotropic  materials.  They  did  this  by  using  the  fact 
a  layered  structure  of  homogeneous  isotropic  material  can  be  treated  as  a  single 
anisotropic  medium  provided  the  layers  are  small  compared  to  wavelength  [40].  For  a 
given  set  of  two  layers  of  material  that  are  sufficiently  thin,  the  effective  permittivities 
are  [100] 
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where  77  =  and  dl  are  the  thicknesses  of  the  two  layers.  Huang  et  al.  used  Equa¬ 
tions  2.34  and  2.35  to  derive  the  appropriate  material  parameters  and  thicknesses  for 
a  cylindrical  cloak  with  parameters  defined  by  Equation  2.33.  They  approximated  the 
radial  variation  in  er  using  ten  anisotropic  layers,  resulting  in  their  isotropic,  homo¬ 
geneous  cloak  having  twenty  layers.  Their  results  showed  good  cloaking  performance 
with  a  12  dB  reduction  in  forward  scattering  when  comparing  a  cloaked  to  an  un¬ 
cloaked  PEC  cylinder.  There  was  a  reduced  scattering  for  all  observation  angles,  but 
the  reduction  was  on  the  order  of  3-4  dB.  Since  Huang’s  method  uses  homogeneous 
materials  in  their  implementation,  the  bandwidth  problem  inherent  with  metamateri¬ 
als  (see  Chapter  III)  is  somewhat  avoided.  However,  extreme  parameter  values  would 
still  be  required  to  manufacture  such  a  cloak,  and  homogeneous  materials  with  these 
extreme  values  over  significant  band  widths  do  not  exist. 


Alii  and  Engheta  have  investigated  an  alternative  method  to  hide  objects  [8,9]. 
Their  method  uses  plasmonic  and  metamaterial  coatings  to  reduce  the  scattering  from 
an  object.  There  are  electrical  size  limitations,  but  the  advantage  of  their  method  is 
the  material  coatings  are  homogeneous  and  isotropic.  A  disadvantage  is  the  plasmonic 
coating  is  dependent  on  the  shape  and  material  properties  of  the  object  to  be  hidden. 


Miller  [65]  proposed  a  method  to  cloak  a  region  of  space  that  uses  sensors  and 
active  surface  sources.  Passive  sensors  are  used  at  the  boundary  of  the  region  to  be 
cloaked.  The  incident  radiation  is  sensed,  and  the  type  of  surface  sources  to  be  placed 
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which  generate  the  appropriate  signals  which  make  the  region  appear  transparent  to 
the  incident  field  to  an  outside  observer  can  be  mathematically  determined.  The 
method  has  the  advantage  of  being  broadband.  The  disadvantage  is  actual  imple¬ 
mentation  for  broad-band  electromagnetic  waves  is  difficult  clue  to  the  speed  of  prop¬ 
agation.  Zambonelli  and  Marnei  proposed  something  similar  by  suggesting  cloaking 
can  be  achieved  by  real-time  scene  manipulation.  Active  sensor  networks  embedded 
in  “garments”  (cloaks)  could  sense  the  incident  electromagnetic  energy  and  then  ra¬ 
diate  the  appropriate  response  to  make  the  cloak  invisible  [105].  This  work  was  a 
theoretical  speculation  as  to  what  might  be  possible  in  the  future  as  active  sensor 
technology  evolves. 

Kildal  et  al.  performed  an  interesting  comparison  between  the  cloaking  method 
described  in  [72]  and  the  efficacy  of  scattering  reduction  when  using  geometric  shaping 
with  hard  and  soft  surfaces.  Specifically,  they  compared  the  theoretical  and  realizable 
performances  of  cloaking  a  PEC  cylinder  of  radius  2 a  using  a  cloak  made  of  metama¬ 
terials  and  a  shaped  geometry  with  a  hard  surface.  The  shaped  geometry  with  hard 
surface  is  limited  by  the  fact  the  incident  angle  of  the  electromagnetic  radiation  must 
be  known,  whereas  a  metamaterial  cloak  theoretically  works  independent  of  angle. 
However,  they  state  this  appears  to  be  the  only  area  where  the  metamaterial  cloak 
is  better  suited,  and  they  are  quick  to  point  out  a  combined  TE/TM  realization  of  a 
shaped  geometry  with  a  hard  surface  is  realizable  and  performs  as  well  as  individual 
TE  or  TM  structures.  Due  to  the  complex  anisotropy  of  the  materials  required  to  im¬ 
plement  the  metamaterial  cloak,  only  a  TM  realization  has  been  physically  realized. 
Additionally  the  metamaterial  cloak  has  an  extremely  narrow  bandwidth,  whereas  the 
shaped  geometry  with  hard  surface  have  approximately  a  20%  operational  bandwidth. 
They  conclude  Pendry  et  al.  ’s  cloaking  method,  while  exotic,  does  not  compare  to 
current  technology  of  shaping  and  RAM  [47]. 
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2. 5  Summary 

This  chapter  presented  background  analysis  showing  cloaking  as  derived  from 
transformation  optics  is  theoretically  possible,  although  bandwidth  will  be  limited  clue 
to  constraints  imposed  by  causality  requirements.  Additionally,  clue  to  singularities 
in  the  required  material  parameters  for  a  two-dimensional  cylindrical  cloak,  and  also 
due  to  the  manufacturing  constraints,  reduced  parameter  sets  were  developed.  The 
reduced  parameter  sets  will  be  further  discussed  in  Chapter  IV.  Alternate  cloaking 
methods  were  discussed,  and  the  work  done  by  Huang  et  al.  will  be  used  in  Chapter 
V. 
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III.  Metamaterials  Background 

The  intent  of  this  chapter  is  not  to  define  or  categorize  metamaterials.  This  is  al¬ 
ready  an  extensive  effort,  and  even  the  mere  definition  of  metamaterial  evokes  much 
discussion  [84],  In  this  effort,  metamaterials  will  be  defined  as  they  were  in  the  work 
by  Lapine  and  Tretyakov  [54]: 

(A)  metamaterial  is  an  arrangement  of  artificial  structural  elements  de¬ 
signed  to  achieve  advantageous  and  unusual  properties. 

No  specific  requirements  regarding  the  material’s  resulting  constitutive  parameters 
are  stated,  nor  does  this  limit  previously  “named”  materials  from  being  included 
in  this  definition.  Examples  of  these  older  materials  are  chiral  materials,  artificial 
dielectrics,  and  artificial  magnetics  [107].  Newer  media  that  have  recently  been  de¬ 
veloped  are  Veselago  media,  which  are  also  known  as  double-negative  (DNG)  media. 
In  a  Veselago  medium,  both  the  effective  permeability  and  permittivity  are  less  than 
zero.  Additionally,  there  are  materials  where  only  the  permittivity  or  permeability 
is  less  than  zero.  These  are  referred  to  as  epsilon-negative  and  mu-negative  materi¬ 
als  respectively  [30].  Note  that  nowhere  in  this  research  are  materials  with  negative 
permeability  or  permittivity  required.  It  must  also  be  stressed  there  is  considerable 
debate  as  to  the  validity  of  the  experimental  results  proving  the  existence  of  DNG 
materials  [66];  however,  the  substance  of  that  debate  it  outside  the  scope  of  this  effort. 

As  stated  in  Section  1.2.2,  metamaterials  are  man-made  materials  with  sub¬ 
wavelength,  often  periodic  structures.  The  structures  are  usually  metallic.  The  ge¬ 
ometry  and  periodicity  of  the  structure  enables  one  to  create  a  material  with  desired 
effective  permittivity  and  permeability  values  that  are  either  isotropic  or  anisotropic. 
This  is  advantageous  because,  as  discussed  in  Chapter  II,  cloaks  require  materials 
with  specific  anisotropic  constitutive  parameters.  Precise  control  of  the  material 
parameters  is  required  during  manufacturing  because  deviations  from  the  specified 
values  result  in  degraded  performance  [42,101].  Therefore,  metamaterials  are  the 
ideal  building  blocks  for  not  only  cloaks,  but  for  any  material  design  created  using 
the  transformation  optics  approach. 
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Metamaterials  are  typically  periodic  lattices  of  unit  cells.  Within  each  unit  cell  is 
some  type  of  metallic  inclusion  whose  shape  determines  the  type  of  bulk  effect  desired. 
In  order  to  realize  bulk  effective  material  parameters,  the  size  of  a  unit  cell,  a  (also 
known  as  the  lattice  constant),  is  much  smaller  than  the  operating  wavelength.  There 
has  not  been  a  theoretically  determined  limit  for  a  such  that  if  the  limit  is  passed, 
the  material  immediately  stops  exhibiting  bulk  constitutive  parameters.  However, 
the  rule  of  thumb  appears  to  be  [86] 


0.01  <  y  <  0.2. 

A 


(3.1) 


This  research  focuses  on  metamaterials  operating  in  the  10  GHz  region.  At  10  GHz, 
the  wavelength  is  three  centimeters.  Therefore,  metamaterial  lattice  constants  will  be 
no  larger  than  six  millimeters.  This  rule  of  thumb  is  crucially  important  if  one  is  to 
properly  use  Maxwell’s  equations  in  their  macroscopic  form  to  analyze  metamaterials. 

A  short-coming  of  metamaterials  is  they  are  dispersive.  This  means  their  con¬ 
stitutive  parameters  vary  as  a  function  of  frequency.  The  clispersiveness  can  be  math¬ 
ematically  represented  using  the  following  general  formulas  [30]: 
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where  tope  is  the  electric  plasma  frequency  and  c opm  is  the  magnetic  plasma  frequency 
which  can  be  controlled  by  varying  the  properties  of  the  unit  cell.  Te  and  Tm  are 
damping  terms  associated  with  losses  in  the  materials. 


What  follows  in  the  subsequent  sections  is  a  discussion  of  some  of  the  more 
common  metamaterial  structures  that  are  used  to  artificially  create  effective  electric 
and  magnetic  mediums.  Note  most  of  the  recent  work  with  metamaterials  has  been 
designing  and  manufacturing  materials  that  will  exhibit  negative  refraction.  As  stated 
above,  negative  refraction  will  not  be  required  for  this  effort,  nor  does  this  effort  intend 
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to  debate  the  existence  of  DNG  materials.  However,  the  structures  used  to  create 
the  desired  constitutive  parameters  will  be  useful.  This  discussion  is  by  no  means 
comprehensive,  and  is  intended  to  give  an  overview  of  how  metamaterials  realize  the 
desired  properties. 

3.1  Creating  Effective  Permittivity 

The  first  work  in  creating  artificial  permittivities  to  result  in  an  electromagnetic 
effect  was  done  by  Kock  in  1946.  He  used  conducting  plates  of  a  certain  shape 
and  prescribed  spacing  to  create  a  microwave  lens  [50].  An  electromagnetic  held 
incident  upon  such  a  structure  will  undergo  a  focusing  effect,  much  the  same  way 
light  does  when  refracted  by  an  optical  dielectric  lens.  In  essence,  the  geometry  of 
the  conducting  plates  created  a  bulk  permittivity  as  seen  by  incident  electromagnetic 
energy.  Similarly,  Kock  used  conducting  paint  on  cellophane  globes  to  create  a  bulk 
permittivity  effect  which  helped  to  reduce  the  weight  of  standard  dielectric  lenses  [51]. 

Since  the  initial  work  by  Kock,  there  have  been  a  great  number  of  contributors 
to  the  held  of  artificial  dielectrics.  Much  effort  was  done  in  the  1950’s  and  1960’s 
to  produce  artihcial  dielectrics  for  use  with  radar.  It  is  now  widely  known  lattice 
structures  of  metallic  spheres,  disks,  or  rods  can  create  a  medium  with  an  effective 
permittivity.  Depending  on  the  spacing  between  the  lattice  objects,  and  the  polariza¬ 
tion  and  frequency  of  the  electromagnetic  held,  the  permittivity  varies  [23].  Rotman 
showed  how  a  periodic  lattice  of  metallic  rods  can  be  used  to  simulate  the  behavior 
of  a  plasma.  Plasmas  have  a  negative  permittivity  while  their  relative  permeability 
is  unity.  The  rodded  artihcial  dielectric  in  [76]  effectively  simulated  a  plasma. 

Pendry  et  al.  used  a  similar  structure  of  thin  wires  to  create  an  effective  negative 
permittivity  in  the  gigahertz  range.  Prior  to  their  work,  the  electric  plasma  frequency 
(Equation  3.2)  was  typically  confined  to  the  optical  frequencies.  They  achieved  the 
frequency  reduction  by  emphasizing  the  requirement  for  very  thin  wires,  which  reduces 
the  electron  density  thereby  lowering  the  plasma  frequency.  Results  were  obtained 
using  numerical  simulations  which  confirmed  the  theoretical  derivations  [71] .  Pendry 
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et  al.  extended  their  results  in  [69]  by  creating  an  array  of  gold-plated  tungsten  wires 
with  a  diameter  of  20  /am.  The  spacing  (lattice  constant)  between  the  wires  was  5 
mm.  As  mentioned  in  Section  2.2,  Gaillot  et  al  used  arrays  of  gold-plated  tungsten 
rods  to  create  the  desired  radial  permittivity  variation  in  a  simulation  of  a  simplified 
cylindrical  cloak  [35]. 

An  alternate  way  to  create  an  artificial  permittivity  is  to  use  an  electric- inductive- 
capacitive  (ELC)  resonator.  An  ELC  resonator  was  created  because  researchers  found 
fabrication  of  a  three  dimensional  wire  array  was  difficult,  and  also  that  extremely 
thin  wires  have  too  high  of  losses  [80] .  An  ELC  resonator  has  a  unit  cell  structure 
shown  in  Figure  3.1. 


Figure  3.1:  ELC  Resonator  [80] 


The  ELC  resonator  operates  as  follows.  For  the  held  configuration  shown  in 
Figure  3.2,  the  current  how  is  due  solely  to  electric  held  coupling.  There  is  no  mag¬ 
netic  held  coupling  because  the  magnetic  held  is  in  the  same  plane  as  the  device. 
Charge  builds  up  on  the  center  T’s  which  creates  a  capacitance.  The  outer  loops  act 
as  two  oppositely  wound  inductors.  This  structure  acts  as  an  equivalent  circuit  where 
a  capacitor  is  in  parallel  with  two  inductors  (Figure  3.2,  lower  right).  For  this  par¬ 
ticular  held  conhguration,  the  capacitive  element  couples  to  the  incident  electric  held 
and  gives  the  metamaterial  the  artihcial  permittivity.  When  the  magnetic  held  is  per¬ 
pendicular  to  the  plane  of  the  device,  as  shown  in  Figure  3.3,  the  currents  in  the  loops 
cancel,  resulting  in  no  magnetic  moment.  There  is  no  coupling  to  the  electric  held  be¬ 
cause  there  is  not  a  significant  capacitance  which  allows  current  to  how.  Controlling 
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d  «  X 


Figure  3.2:  ELC  Electric  Field  Coupling 


Figure  3.3:  ELC  Operation  with  No  Magnetic  Field  Coupling 

the  capacitance  and  the  inductance  is  done  by  varying  the  parameters  a,  d,  l ,  w,  and 
g  shown  in  Figure  3.1.  These  parameters  determine  the  resonant  frequency  of  the 
device.  Schurig  et  al.  obtained  simulated  S-parameter  measurements  on  an  ELC  ar¬ 
ray  using  Microwave  Design  Studio,  an  electromagnetics  simulation  software  package 
based  on  a  finite  integral  time  domain  formulation.  They  also  built  a  single-layered 
array  of  structures  and  obtained  excellent  agreement  between  extracted  parameters 
(see  Section  3.3)  from  simulated  and  experimentally  obtained  measurements. 
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3.2  Creating  Effective  Permeability 

There  was  a  significant  push  in  the  late  1990’s  to  create  magnetic  effects  from 
non-magnetic  materials  due  to  the  quest  to  create  DNG  media.  A  variety  of  structures 
has  been  analyzed  which  generate  magnetic  responses.  The  following  subsections  dis¬ 
cuss  these  structures.  The  primary  focus  of  the  research  has  been  manufacturing 
a  DNG  medium.  However,  these  various  structures  are  beneficial  because  they  en¬ 
able  the  ability  to  create  metamaterials  with  a  specified  permeability  at  a  desired 
frequency.  The  different  designs  each  have  their  own  positive  and  negative  aspects. 

3.2.1  Edge  Coupled  Split  Ring  Resonator.  Pendry  et  al.  developed  what  is 
commonly  called  the  edge  coupled  split  ring  resonator  (EC-SRR).  The  basic  building 
block  for  the  EC-SRR  is  the  split  ring  structure  shown  in  Figure  3.4.  The  gold 
structures  in  the  figure  are  thin  layers  of  metal  mounted  on  a  substrate.  The  device 
is  excited  by  a  plane  wave  propagating  in  the  plane  of  the  device  with  the  magnetic 
field  normal  to  the  plane  of  the  rings.  Since  the  magnetic  flux  is  perpendicular  to  the 


Figure  3.4:  Edge  Coupled  Split  Ring  Resonator 
plane  of  the  rings,  an  electromotive  force  results  causing  current  to  flow.  In  general, 


39 


the  electromotive  force  is  [74] 


EMF  =  j)E-dl  =  -ju  J  B  ■  dS.  (3.4) 

The  average  radius  of  the  two  rings  is  rg.  Thus,  EMF  =  — jumr^Bz  [62],  The 
flowing  current  in  the  inner  and  outer  rings  creates  the  magnetic  moment,  resulting 
in  an  artificial  permeability.  Current  flow  is  possible  dne  to  the  capacitance  in  the 
ring  gaps  (t  and  d  in  Figure  3.4).  A  stronger  capacitance  will  result  in  more  current 
flowing  because  capacitance,  charge,  and  voltage  are  related  by  [74] 

C=%.  (3.5) 

Current  flowing  in  the  outer  and  inner  loops  results  in  the  device  having  an  inductance 
which  is  related  to  the  radius  of  the  rings.  The  capacitance  and  inductance  results  in 
the  device  having  a  resonant  frequency  of  -7=.  Detailed  analytic  derivations  for  C 
and  L  are  provided  in  [78]. 

Typically,  an  EC-SRR  is  called  an  SRR.  Such  a  configuration  of  an  EC-SRR 
results  in  a  material  with  an  effective  permeability  described  by  Equation  3.3,  where, 
for  this  geometry,  topm  and  Ym  are  controlled  by  the  various  dimensions  of  the  geom¬ 
etry  within  the  unit  cell  in  addition  to  the  lattice  spacing  [70].  Aydin  et  al.  showed 
using  simulations  and  experimentally  that  increasing  both  d  and  t  shifts  the  resonant 
frequency  higher  due  to  the  reduced  capacitance  [12].  An  array  of  SRRs  as  shown  in 
Figure  3.4  in  combination  with  an  array  of  metallic  wires  was  used  to  create  the  first 
composite  DNG  medium  [88]  (although  as  stated  earlier,  there  is  some  debating  the 
meaning  of  these  results). 

3.2.2  Omega  Split  Ring  Resonator.  Another  geometry  used  to  create  mag¬ 
netic  effects  is  the  f!  ring.  However,  unlike  the  EC-SRR,  which  has  no  electric  effects 
without  being  used  in  conjunction  with  a  rod  lattice,  the  D  ring  couples  to  the  electric 
and  magnetic  fields.  A  typical  f!  ring  geometry  is  shown  in  Figure  3.5.  Currents  are 
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Figure  3.5:  fi-Riug  Geometry 

excited  in  the  rod-portion  of  the  hi  ring  due  to  the  incident  electric  held.  Charge  builds 
up  in  the  gap  and  the  rod  ends,  resulting  electric  moments  which  creates  the  artificial 
permittivity.  The  magnetic  held  is  normal  to  the  plane  of  the  device.  This  excites 
a  current,  which  hows  in  the  ring  portion  of  the  geometry,  creating  the  magnetic 
moments  resulting  in  an  artihcial  permeability.  An  array  of  Q  rings  is  bianisotropic. 
This  is  because  the  charge  build  up  in  the  rings  resulting  from  the  current  induced 
by  the  magnetic  held  results  in  a  series  of  electric  dipoles.  Simovski  and  ffe  showed 
that  two  hi  geometries  printed  in  opposite  directions  on  opposing  sides  of  a  dielectric 
do  not  have  any  magnetoelectric  [85].  The  opposing  direction  of  current  how  in  the 
rings  results  in  the  electric  dipoles  induced  by  the  magnetic  held  cancelling.  Note  the 
original  EC-SRR  developed  by  Pendry  et  al.  was  also  shown  to  have  bianisotropic 
effects  [62],  although  these  effects  can  be  minimized  by  varying  the  EC-SRR  orienta¬ 
tion  throughout  a  lattice.  Huangfu  et  al.  used  an  O-like  structure  to  build  a  DNG 
medium  whose  effects  manifested  from  12-13  GHz.  Only  a  periodic  structure  consist¬ 
ing  of  elements  shown  in  Figure  3.5  were  required.  They  did  not  have  to  include  a 
lattice  of  thin  metallic  wires  to  generate  the  negative  permittivity.  This  is  because  the 
rods  connecting  the  fl  structures  serve  this  function.  Huangfu  et  aV s  results  showed 
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more  bandwidth  where  the  DNG  effect  took  place  [41],  and  they  also  reported  smaller 
losses  than  those  experiments  which  used  ring  geometries. 

3.2.3  Additional  Structures.  O’Brien  and  Pendry  developed  what  is  now 
called  an  axially  symmetric  ring.  The  geometry  is  very  similar  to  the  EC-SRR  in  the 
previous  subsection  except  instead  of  rings,  rectangular  structures  are  used.  This  is 
shown  in  Figure  3.6.  The  advantage  of  this  type  of  geometry  is  the  structure  is  more 


Figure  3.6:  Axially  Symmetric  Ring 

conducive  to  manufacturing  at  the  sub-micron  wavelengths  [68].  These  smaller  struc¬ 
tures  would  enable  the  desired  effects  to  manifest  at  optical  frequencies  provided  the 
metal  maintains  its  good  conductor  properties  at  these  higher  frequencies.  O’Brien 
and  Pendry  used  a  such  a  structure  in  conjunction  with  a  lattice  of  thin  metallic  wires 
to  create  a  DNG  medium  operational  at  76  THz.  There  are  also  results  presented 
in  [30]  where  a  negative  permeability  was  demonstrated  using  an  array  of  axially 
symmetric  SRRs  between  8.2  and  8.7  GHz. 

An  additional  metamaterial  structure  is  called  the  S-ring.  Rs  geometry  is  shown 
in  Figure  3.7.  Much  like  the  Q-ring,  the  S-ring  geometry  does  not  require  the  use  of 
metallic  rods  to  generate  an  electric  effect.  Chen  et  al.  designed  and  fabricated  a  DNG 
metamaterial  using  the  S-ring  structures,  with  an  operational  frequency  of  10.9  -  13.5 
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GHz,  which  is  a  significant  increase  in  bandwidth  compared  to  early  implementations 
[17].  They  performed  a  detailed  theoretical  analysis  and  develop  the  theory  which 
governs  the  effective  constitutive  parameters  of  an  S-ring  metamaterial  [16]. 

There  are  other  types  of  structures  that  can  be  used.  In  fact,  when  Schurig  et  al. 
manufactured  the  first  simplified  cylindrical  cloak,  the  unit  cells  contained  aspects  of 
the  axially  symmetric  SRR  in  addition  to  the  ELC-resonator.  A  picture  of  the  basic 
geometry  of  the  unit  cell  is  shown  in  Figure  3.8.  Changing  various  dimensions  of 
the  cell  (7,  s,  w,  and  r)  alters  the  cell  capacitance  and  inductance,  thereby  creating 
the  desired  permittivity  and  permeability  specified  by  the  simplified  cloak’s  material 
parameters. 

Obviously,  there  are  limitless  types  of  geometries  which  could  be  used  to  create 
artificial  electric  and  magnetic  effects.  However,  in  order  avoid  simply  trial  and  error, 
one  must  understand  how  the  fields  interact  with  the  devices. 

3.3  Measuring  Metamaterial  Constitutive  Parameters 

As  stated  in  the  previous  sections,  there  has  been  much  work  developing  peri¬ 
odic  lattices  of  unit  cells  which  exhibit  artificial  permittivity  and  permeability.  This 
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Figure  3.8:  Unit  Cell  Geometry  for  Simplified  Cloak  Construction  [79] 

has  led  to  a  good  understanding  of  how  to  create  bulk  artificial  electromagnetic  ef¬ 
fects  using  metamaterials.  For  most  experiments,  the  proof  the  metamaterial  was 
exhibiting  DNG  behavior  was  obtained  by  empirically  measuring  the  refracted  angle 
of  incident  energy  at  various  frequencies.  If  this  angle  was  negative,  it  was  concluded 
the  material  had  a  negative  index  of  refraction  and  therefore  simultaneously  negative 
permittivity  and  permeability.  This  can  be  considered  an  indirect  measurement  of 
the  constitutive  parameters  because  the  materials  themselves  were  not  explicitly  mea¬ 
sured.  Therefore,  direct  methods  to  measure  the  permittivity  and  permeability  have 
been  developed.  There  has  been  much  discussion  during  the  past  decade  whether 
standard  material  parameter  retrieval  algorithms  can  be  used  when  characterizing 
metamaterials.  In  the  following  sections,  a  standard  material  retrieval  algorithm  is 
presented  followed  by  various  analyses  done  on  metamaterial  characterization  using 
this  algorithm.  Limitations  of  the  retrieval  algorithm  are  discussed  with  possible 
alternative  measurement  techniques  documented. 

3. 3. 1  Nicolson-Ross-  Weir  Algorithm.  A  common  method  to  extract  permit¬ 
tivity  and  permeability  from  a  homogeneous  material  slab  is  to  obtain  S-parameter 
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measurements  and  then  use  the  Nicolson-Ross-Weir  (NRW)  algorithm  to  extract  the 
constitutive  parameters  [67,99].  The  NRW  algorithm  is  derived  as  follows.  From  a 
sample  of  an  infinite  length  of  planar  material,  it  can  be  shown  that 


R[l  —  P 2 
1  -  R2P2  ’ 


(3.6) 


P  (1  -  II- 
1  -  Wp 5  ’ 


(3.7) 


where  P  =  e~^kd  and  is  called  the  phase  delay,  k  is  the  wave  number  in  the  slab  of 
material,  and  d  is  the  material  slab  thickness  [39].  The  reflection  coefficient,  R ,  from 
an  infinite  slab  of  planar  material  in  free  space  is 


(3.8) 


where  Z  and  ZQ  are  the  intrinsic  impedance  of  the  material  and  free  space  respec¬ 
tively  and  Z  =  y/jtje  [13].  It  has  been  assumed  a  transverse  electromagnetic  wave 
propagating  in  its  fundamental  mode  will  be  used  to  interrogate  the  material.  Note 
also  that  k  =  cuy fjle  =  u>y/n0e 0/+er  =  -y' /irer  where  c  is  the  free  space  speed  of  light. 
Thus,  by  finding  R  and  P ,  one  can  use  their  respective  definitions  to  find  /ir  and  er. 

First,  one  must  find  expressions  for  R  and  P  using  only  the  S-parameters.  It  is 
possible  to  solve  Equations  3.6  and  3.7  for  P  [39]. 


P 


S' 
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1  -  RS n 


(3.9) 


Equation  3.9  can  be  used  to  write  a  quadratic  expression  for  R. 


R2  -  2 QR  +  1  =  0 


(3.10) 


where 


Q  = 


si  i 


5|i  +  1 


2Sh 


(3.11) 
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Therefore  R  is  found  to  have  two  solutions. 


Ri.2  =  q±  VQ 2  -  1, 


(3.12) 


The  sign  choice  is  made  such  that  \R\  <  1  [67].  Note  R  is  now  expressed  simply 
in  terms  of  the  S n  and  S2 1  measurements.  Once  R  is  found,  P  can  be  found  using 
Equation  3.9  [39]. 

Equation  3.8  can  be  expanded  using  the  definition  of  Z . 


R 


£r£o 


£r&o 


+ 


This  can  be  simplified  to 


Rearranging  Equation  3.14  leads  to  the  final  form 


z  = 


1  +  R 
1  -R' 


(3.13) 


(3.14) 


(3.15) 


where  z  is  commonly  called  the  normalized  intrinsic  impedance.  Similarly,  the  ex¬ 
panded  expression  for  k  can  be  used  to  expand  the  expression  for  P  as 


P  =  e~jkd  = 


(3.16) 


Solving  Equation  3.16  for  yjnrer  yields  [67] 

, _  iclnP 

y  =  ^J[irer  = - .  (3.17) 

(jJ 

Using  Equations  3.15  and  3.17,  one  can  write  expressions  for  /ir  and  er  in  terms  of  R 
and  P  which  can  be  expressed  using  only  the  measured  S-parameters  (Equations  3.9 
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and  3.12)  [67], 


Mr  =  yz, 


(3.18) 


y 

z 


Thus,  nr  and  er  can  be  found  using  only  the  S-parameter  measurements. 


3.3.2  Metamaterial  Constitutive  Parameter  Extraction.  Smith  et  al.  were 
the  first  to  use  the  NRW  algorithm  to  extract  constitutive  parameters  from  metamate¬ 
rial  measurements  obtained  from  simulations.  They  used  reflection  and  transmission 
coefficients  (effectively  S\  i  and  S21)  generated  using  a  transfer  matrix  simulation  on  a 
periodic  array  of  wires,  a  periodic  array  of  SRRs,  and  a  periodic  array  of  a  wire-SRR 
structure.  The  extracted  constitutive  parameters  for  these  media  were  consistent 
with  Equations  3.2  and  3.3  [89].  Markos  and  Soukoulis  obtained  similar  results  us¬ 
ing  simulated  reflection  and  transmission  coefficient  data  to  obtain  the  constitutive 
parameters  of  metamaterials  [61]. 

There  are  some  issues  which  must  be  considered  when  extracting  constitutive 
parameters  from  metamaterials  using  the  NRW  algorithm.  Smith  et  al.  emphasize 
the  process  is  only  accurate  when  measuring  metamaterials  which  do  not  have  chiral 
or  bianisotropic  effects  [89].  A  metamaterial  consisting  of  an  array  of  wires  does 
not  exhibit  any  of  these  behaviors.  However,  as  shown  by  Marques  et  al. ,  an  array 
of  SRRs  does  exhibit  bianisotropic  effects  [62],  This  does  not  mean  Smith  et  al.' s 
material  parameters  are  not  accurate  because  their  devices  possessed  small  chiral 
effects.  Thus,  admittedly  their  measurements  could  not  be  considered  exact,  but 
they  believed  they  were  still  a  good  estimate  of  the  constitutive  parameters  [89]. 
Thus,  when  performing  these  types  of  measurements,  one  must  have  a  general  idea 
as  to  the  extent  the  material  will  exhibit  bianisotropic  behavior.  If  it  is  theoretically 
believed  the  material  will  have  large  chiral  effects,  than  an  alternate  method  to  extract 
material  parameters  will  be  required  [58].  This  is  due  to  the  fact  the  wave  number, 
k,  is  a  function  of  not  only  /x  and  £  but  also  the  magnetoelectric  coupling  coefficients. 
Hence,  the  NRW  algorithm  assumes  the  material  being  analyzed  is  isotropic  and 
homogeneous. 
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Another  issue  encountered  when  extracting  material  parameters  from  S-parameter 
measurements  of  metamaterials  using  the  NRW  algorithm  has  to  do  with  the  sign 
choice  on  the  square  root  used  in  the  parameter  extraction.  To  help  explain,  consider 
an  alternate  form  of  the  equations  used  to  determine  z.  As  noted  in  the  previous  sec¬ 
tion,  z  is  the  normalized  impedance.  Equation  3.15  can  be  expanded  and  expressed 
using  only  the  S-parameter  measurements. 


(1  -  Sn)2  - 


(3.19) 


Note  that  z  will  have  both  a  real  and  imaginary  component  due  to  the  complex 
nature  of  the  S-parameter  measurements.  For  passive  materials,  the  real  part  of  the 
normalized  impedance  must  be  positive  [22],  which  makes  the  sign  choice  in  Equation 
3.19  obvious. 

The  assumption  of  passive  materials  also  helps  in  determining  another  sign 
choice  in  the  parameter  extraction.  Recall  that  P  =  and  that  k  =  oj^iirix0£r£0. 
The  index  of  refraction,  n  can  be  defined  as 


n  =  nrer.  (3.20) 

Additionally,  the  free  space  wave  number,  ka  is  defined  as 


K  ^  \J  1-1,-Po-  (3.21) 

Using  Equations  3.20  and  3.21,  the  equation  defining  the  phase  advance,  P,  can  be 
rewritten  as 

P  =  e-jk0nd .  (3.22) 

Equation  3.22  can  be  expressed  solely  in  terms  of  the  S-parameter  measurements. 

e-jk0nd  =  x  ±  _  x2  (3.23) 
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where 


1 


X 


2s2l{i-si1  +  si1y 


(3.24) 


The  sign  choice  in  Equation  3.22  is  made  knowing  the  imaginary  part  of  the  refractive 
index  for  passive  materials  must  be  less  than  zero  [89].  Solving  for  the  real  part  of  the 
refractive  index  does  present  problems  clue  to  the  multiple  branches  of  the  complex 
logarithm.  Specifically,  when  solving  Equation  3.22,  one  gets 


In  P  +  jm2n  =  — jk0nd , 


(3.25) 


where  the  jm2n  term  accounts  for  the  branch  choice.  Note  that  In  P  will  have  both 
a  real  and  imaginary  component  because  P  is  a  complex  quantity.  Hence,  this  can 
be  rewritten  as 

Re  [ln(P)]  +  jlrn  [ln(P)]  +  jrn2n  =  —jk0nd.  (3.26) 

Solving  for  n  yields 

n  =  y —  ([— Im  ln(P)]  —  2rmr  +  jRe  [ln(P)]) .  (3.27) 

Due  to  branch  cuts,  the  real  part  of  n  is  ambiguous.  This  ambiguity  typically  does 
not  cause  problems  because  the  thickness  of  a  material  sample  is  usually  known.  This 
may  not  be  the  case  for  a  metamaterial.  Ambiguities  can  arise  because  the  sample 
width  and  the  reference  plane  (first  effective  boundary)  of  a  sample  of  metamaterial 
are  sometimes  not  known.  The  reference  plane  is  the  location  after  which  reflected 
waves  from  plane  wave  incidence  exhibit  plane  wave  behavior  [22],  Additionally,  the 
sample  width  of  a  metamaterial  may  also  be  ambiguous  clue  to  the  fact  metamate¬ 
rials  may  not  have  well  dehned  surfaces.  In  [89],  Smith  et  al.  note  that  more  than 
one  sample  thickness  must  be  measured  in  order  to  identify  the  correct  branches. 
Different  thicknesses  should  result  in  the  same  material  parameters;  thus,  the  two 
measurements  are  compared  and  the  branch  which  makes  the  solutions  identical  is 
chosen.  The  authors  suggest  using  sample  thicknesses  as  small  as  possible  in  order 
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to  make  the  branch  cuts  further  apart,  which  will  make  the  selection  of  the  correct 
branch  easier. 

Chen  et  al.  present  an  improved  method  for  choosing  the  correct  branch  cut. 
They  use  the  requirement  the  constitutive  material  parameters  are  continuous  func¬ 
tions  of  frequency  [89]  and  use  an  iterative  approach.  They  assume  they  have  a  correct 
value  for  the  refractive  index  at  a  given  frequency  and  use  a  Taylor  series  to  expand 
P  in  Equation  3.22  at  the  next  sampled  frequency  and  choose  the  branch  cut  which 
as  close  to  possible  enforces  the  continuity  requirement  [22], 

Other  authors  have  used  the  NRW  algorithm  to  extract  constitutive  material 
parameters  from  both  simulated  and  experimental  S-parameter  data  from  metama¬ 
terials.  Ziolkowski  simulated  arrays  of  SRRs  and  also  arrays  of  capacitively  loaded 
strips.  His  designs  were  simulated  using  High  Frequency  Structure  Simulator  (HFSS), 
a  commercial  FEM-based  electromagnetics  software  package  and  also  using  the  com¬ 
mercial  FDTD  package  produced  by  Ocotillo  ElectroMagnetics.  He  provided  a  slight 
deviation  to  the  extraction  formulas  because  he  noted  when  SR  and  SR  are  nearly 
zero,  choosing  the  sign  of  a  square  root  becomes  ambiguous.  His  measured  results 
compared  favorably  to  experimental  results,  with  the  errors  attributed  to  imprecision 
in  the  manufacturing  of  the  devices  [108]. 

Greegor  et  al.  simulated  a  standard  SRR  and  wire  configuration  in  the  13  -  15 
GHz  range  using  Microwave  Design  Studio  .  They  used  the  simulated  S-parameters 
to  determine  the  refractive  index,  n.  They  then  manufactured  the  metamaterial 
and  obtained  measurements.  Values  for  n  obtained  from  the  simulated  data  were 
in  agreement  within  20%  of  the  values  obtained  using  measured  data  [37].  Smith 
et  al.  used  the  same  parameter  retrieval  method  they  developed  in  [89]  to  extract 
the  material  parameters  from  metamaterials  whose  unit  cells  were  not  homogeneous. 
They  simulated  metal  wire  and  SRR  unit  cell  arrays  using  HFSS.  They  found  no 
changes  were  necessary  provided  the  unit  cells  are  periodic  along  the  direction  of 
propagation  [90]. 
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Chen  et  al.  used  the  same  NRW  algorithm,  although  their  equations  were 
slightly  different  because  they  interrogated  the  material  in  a  waveguide  using  a  dom¬ 
inant  TEio  mode,  which  changes  the  impedance  of  the  medium  and  the  resulting 
equations  for  /ir  and  er  developed  in  Section  3.3.1.  They  obtained  S-parameter  mea¬ 
surements  of  an  array  on  metamaterial  unit  cells  with  an  SRR-type  geometry.  By 
rotating  the  geometry  of  the  SRRs  and  then  performing  different  S-parameter  mea¬ 
surements,  Chen  et  al.  were  able  to  extract  the  entire  permittivity  and  permeability 
tensors,  something  which  had  not  yet  been  done  using  measured  data  [19]. 

For  the  work  discussed  in  Chapter  VI,  the  parameter  retrieval  method  discussed 
in  this  and  the  previous  subsections  are  used.  The  software  package  to  perform  simu¬ 
lated  S-parameter  measurements  on  unit  cells  used  was  Cornsol  Multiphysics,  which 
has  not  been  cited  as  used  for  this  type  of  work  in  the  literature.  However,  at  the  2008 
Cornsol  Multiphysics  conference,  Urzhumov  presented  a  detailed  procedure  to  extract 
constitutive  parameters  from  simulated  S-parameter  measurements  using  Cornsol.  His 
method  mirrored  that  presented  in  this  subsection,  with  the  only  differences  being  in 
the  implementation  in  the  software  [92] .  This  will  be  discussed  further  in  Chapter  VI 

3-4  Alternate  Parameter  Retrieval  Method 

An  alternate  technique  to  extract  material  parameters  for  metamaterial  unit 
cells  was  put  forth  by  Smith  and  Pendry.  The  process  is  called  field  averaging.  The 
extraction  of  material  parameters  using  this  method  can  only  be  accomplished  from 
simulations  of  unit  cells.  It  is  not  applicable  to  experimental  data.  Fundamentally, 
field  averaging  is  a  rather  simple  concept.  It  uses  the  basic  relationship  between  the 
electric  (magnetic)  field  and  electric  (magnetic)  displacement  vectors.  Specifically, 
for  materials  with  no  magnetoelectric  coupling,  the  relationships  are 

D  =  e0£r  ■  E  B  =  yU0/Tr  •  H  (3.28) 
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A  plane  wave  excitation  in  one  of  the  unit  axis  directions  is  assumed.  Additionally,  the 
medium  is  assumed  to  be  reciprocal  i.e.  the  constitutive  parameter  tensors  only  have 
nonzero  diagonal  terms.  It  is  possible  to  solve  for  the  component  values  ( ex,£y,£z , 
Hx,  ny,  Hz)  by  using  the  simulated  values  for  all  held  components  and  performing  a 
simple  division.  As  an  example,  for  a  plane  wave  excitation  that  has  only  an  ^-directed 
electric  held,  ex  = 

In  a  metamaterial  unit  cell  simulation,  the  held  values  are  known  at  all  points 
in  the  geometric  grid  used  to  model  the  material.  Obviously  one  cannot  simply 
perform  this  division  at  each  location  in  the  grid  because  the  goal  is  to  define  bulk 
material  parameters.  Thus,  Smith  and  Pendry  define  held  averages  along  the  various 
faces  and  edges  of  the  unit  cells,  which  are  then  used  as  the  bulk  values  to  calculate 
the  permittivity  and  permeability  tensors  [91].  Simulations  of  various  metamaterial 
cell  geometries  were  performed  and  material  parameters  extracted.  Results  closely 
matched  those  obtained  using  the  NRW  retrieval  techniques  previously  discussed.  It 
was  noted  during  the  analysis  the  extracted  parameters  exhibited  a  spatial  dispersion. 
If  the  unit  cells  were  filled  with  nothing  i.e.  a  metamaterial  composed  solely  of  free 
space,  the  extracted  material  parameters  were  /i0  and  eQ  but  multiplied  by  a  sm^  - 
type  term.  Thus,  the  spatial  dispersion  in  the  parameters  was  determined  to  be  a 
function  of  the  simulation.  Fortunately,  it  appears  the  dispersion  can  be  eliminated 
simply  by  scaling  by  the  term  [91].  The  advantage  of  parameter  extraction  using 
held  averaging  is  there  are  no  decisions  required  in  terms  of  signs  on  square  roots 
or  on  logarithm  branches.  The  disadvantage  is,  although  quite  simple  in  concept, 
implementation  is  fairly  complex  and  computationally  expensive.  There  are  other 
approaches  similar  to  held  averaging,  such  as  held  summation  that  are  conceptually 
similar,  the  exception  being  how  the  bulk  held  values  are  computed  [57]. 

3.5  Tunable  Metamaterials 

As  discussed  in  Sections  3.1  and  3.2,  it  is  possible  to  create  metamaterials  with 
specihc  electric  and  magnetic  properties.  However,  once  the  devices  are  created,  the 
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resonant  frequencies  are  fixed.  There  has  been  research  into  unique  ways  to  alter  the 
resonant  frequency  of  general  SRR-structures  following  fabrication  of  the  devices. 

Aydin  and  Ozbay  manufactured  SRRs,  similar  to  the  one  showed  in  Figure  3.4, 
and  altered  the  SRR’s  capacitance  by  mounting  capacitors  between  various  aspects 
of  the  SRR.  This  is  shown  in  Figure  3.9.  Note  how  the  capacitors  are  mounted  both 


Figure  3.9:  SRR  devices  with  capacitors  mounted  post-fabrication  [11] 

between  the  gaps  in  each  of  the  rings  and  also  between  the  inner  and  outer  rings 
themselves.  The  mounted  capacitor  obviously  changed  the  overall  capacitance  of  the 
cell,  which  would  down  shift  the  magnetic  plasma  frequency  ( ujprn  in  Equation  3.3). 
Aydin  and  Ozbay  proved  the  change  in  frequency  by  measuring  the  S-parameters 
of  the  device  before  and  after  capacitor  mounting  [11],  As  expected,  when  different 
capacitor  values  were  used,  different  magnitudes  of  frequency  shift  were  observed. 

Shadrivov  et  al.  did  similar  work,  except  instead  of  mounting  a  standard  capac¬ 
itor  between  the  gaps  or  rings  of  a  standard  SRR,  they  mounted  a  variable-capacitor- 
cliode,  whose  capacitance  can  be  controlled  by  applying  a  DC  bias  voltage.  Their 
device  is  shown  in  Figure  3.10.  By  varying  the  voltage  from  -1  to  1  volt,  they  were 
able  to  achieve  a  tuning  range  of  630  MHz  [82],  Similar  work  was  done  by  Gil  et  al. 
who  mounted  a  varactor  onto  a  SRR  and  achieved  approximately  500  MHz  of  tuning 
capability  [36] .  Hand  and  Cummer  also  tuned  the  resonant  frequency  of  an  SRR,  but 
they  used  a  microelectromechanical  systems  (MEMS)  switch  to  alter  the  capacitance 
of  the  structure.  Unlike  Gil  et  al.  and  Shadrivov  et  al,  Hand  and  Cummer  did  not 
obtain  a  tuned  range  of  operating  frequencies.  Rather,  turning  the  switch  on  or  off 
created  two  different  resonant  frequencies  approximately  800  MHz  apart  [38]. 
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Figure  3.10:  Variable-capacitor-diode  mounted  on  an  SRR  device  [82] 

Degiron  et  al.  performed  conceptually  similar  experiments  in  that  their  desire 
was  to  control  the  magnetic  resonance  of  an  SRR-type  structure.  However,  their 
approach  was  rather  unique.  They  used  a  single  gapped  ring  of  copper  (in  essence, 
only  one  of  the  rings  of  an  edge-coupled  SRR)  and  placed  a  piece  of  n-type  silicon  in 
the  gap.  This  is  shown  in  Figure  3.11.  An  815  nrn  laser  was  used  to  illuminate  the 


Figure  3.11:  N-type  silicon  in  SRR  gap 


silicon  at  intensity  levels  from  0-5  mW.  S-parameter  measurements  were  taken,  with 
the  primary  focus  being  on  the  S^i  he.  transmission  parameter.  As  the  intensity  of 
the  illumination  increased,  the  conductivity  of  the  silicon  increased,  thereby  shorting 
out  the  ring  and  removing  the  resonance,  which  manifested  itself  as  an  increase  in  the 
5'2i  measurement.  Degiron  et  al.  also  found  controlling  the  gap  between  the  silicon 
and  the  edge  of  the  rings  had  an  effect  on  the  resonant  frequency  of  the  device  [29]. 
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Finally,  Rederus  used  a  variable  capacitive  MEMS  device  mounted  in  gap  in  an 
SRR.  The  MEMS  device  consisted  of  cantilevered  beams  that,  when  activated,  would 
alter  the  capacitance  of  the  SRR,  thus  changing  the  resonant  frequency  [75].  The 
main  advantage  of  Rederus’s  design  was  his  device  was  fabricated  in-situ  with  the 
SRR,  resulting  in  a  much  cleaner  device  than  those  fabricated  by  Aydin  and  Ozbay 
and  Shadrivov  et  al.  Due  to  technical  issues,  S-parameter  measurements  were  not 
taken,  but  it  was  shown  the  activated  MEMS  device  could  add  0.54  -  0.62  pF  of 
capacitance  to  the  SRR. 

3. 6  Summary 

This  chapter  presented  various  unit  cell  structures  which  have  been  used  to 
create  effective  bulk  permittivity  and/or  permeability  effects  in  man-made  materials. 
The  NRW  algorithm,  which  can  be  used  to  extract  the  material  parameters  from  some 
geometries,  was  documented  with  its  limitations  and  ambiguities  when  measuring 
metamaterials  were  delineated. 
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IV.  Improved  Simplified  Parameters  for  Two-Dimensional 

Cylindrical  Cloaks 

In  Section  2.2,  it  was  shown  Schurig  et  al.  made  an  error  when  deriving  the  wave 
equation  for  a  TMZ  electromagnetic  field  in  an  anisotropic  cloak.  Their  error  was 
assuming  fig  was  constant  [101].  However,  as  clearly  seen  in  Equation  1.21,  the  ideal 
value  for  fig  is  not  spatially  invariant.  Thus,  the  simplified  parameter  sets  developed, 
while  shown  to  work  reasonably  well,  were  not  correctly  obtained.  In  this  chapter, 
all  constraint  equations  for  the  material  parameters  of  an  ideal  cylindrical  cloak  are 
derived  using  the  correct  form  of  the  wave  equation.  Using  these  constraint  equations, 
it  is  first  shown  how  varying  fig  can  control  the  amount  of  field  transmitted  into  the 
cloak’s  hidden  region.  It  is  then  shown  how  all  constraint  equations  can  be  used  to 
derive  simplified  material  parameter  sets  whose  complexity  can  be  tailored  depending 
on  the  available  manufacturing  capabilities  of  metamaterials. 


4-1  Constraint  Equations 

In  order  to  develop  all  the  constraint  equations  on  the  material  properties  of 
an  ideal  cylindrical  cloak,  the  correct  wave  equation  must  first  be  derived.  Since  it  is 
known  a  priori  the  ideal  parameters  are  ^-invariant  but  not  r-invariant,  the  general 
wave  equation  for  TM2  fields  in  anisotropic  media  shown  in  Equation  2.6  can  be 
expanded  analytically  to 
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Recall  the  ideal  parameters  for  an  ideal  cylindrical  cloak  for  TM2  incident  fields 


are 
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When  one  uses  the  ideal  cylindrical  cloak’s  material  parameters  shown  in  Equation  4.3 
in  the  general  wave  equation  shown  in  Equation  4.1,  the  result  is  the  wave  equation 
for  TM2  fields  in  an  ideal  cylindrical  cloak. 
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As  discussed  in  Section  2.2,  the  original  simplified  material  parameters  for  cylindrical 
cloaks  for  TMZ  incident  field  are 
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The  material  parameters  shown  in  Equation  4.5  were  initially  thought  to  satisfy  the 
same  wave  equation  as  the  ideal  parameter  set.  However,  Yan  et  al.  showed  the 
procedure  leading  to  this  conclusion  was  questionable  by  proving  the  simplified  and 
ideal  parameter  sets  satisfy  different  wave  equations  [101].  This  is  explicitly  seen  by 
substituting  the  simplified  material  parameters  shown  in  Equation  4.5  into  Equation 
4.1.  The  resulting  wave  equation  is  shown  below. 
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Note  the  subtle  difference  in  Equation  4.6  compared  to  Equation  4.4.  The  wave 
equation  using  the  original  simplified  parameters  has  a  1  jr  factor  in  front  of  the 
dEz/dr  term.  The  wave  equation  using  the  ideal  parameters  has  a  factor  of  l/(r  —  a ) 
for  this  same  term.  Thus,  for  values  such  that  r  3>  q,  the  field  behavior  of  the  two 
cloaks  will  be  similar  [101],  but  they  certainly  do  not  satisfy  the  same  wave  equation. 
Since  the  ideal  cylindrical  cloak  effectively  guides  electromagnetic  energy  around  the 
hidden  region,  it  makes  sense  that  a  simplified  cloak,  whose  wave  equation  differs  from 
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that  of  the  ideal  cloak,  has  energy  transmitted  into  this  same  region  [101].  It  also 
makes  sense  that,  because  the  wave  equations  are  similar  for  the  region  away  from 
r  =  a,  the  simplified  parameter  set  performed  reasonably  well  in  terms  of  cloaking 
capability. 


By  comparing  Equations  4.1  and  4.4,  one  finds  the  cylindrical  cloak’s  constitu¬ 
tive  parameters  must  meet  the  following  conditions  in  order  to  achieve  ideal  cloaking 
for  TMZ  fields  [63], 
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The  third  constraint  equation  given  in  Equation  4.9  had  not  appeared  in  the  literature 
prior  to  this  research  and  forms  the  basis  for  the  alternative  simplified  parameters 
proposed  later  in  this  chapter. 


The  ideal  cylindrical  cloak  also  has  an  impedance  at  the  outer  boundary,  r  =  b, 
that  matches  free  space. 
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The  simplified  parameters  shown  in  Equation  4.5  satisfy  Equations  4.7  and  4.8 
but  do  not  satisfy  Equation  4.9.  Additionally,  the  impedance  mismatch  at  r  —  b 
for  the  cylindrical  cloak  with  the  same  simplified  constitutive  parameters  and  using 
a  —  A  and  b  =  2A  is 


=  0.5  (4.11) 

r=b 

Obviously,  there  will  be  a  scattered  field  from  the  simplified  cylindrical  cloak  with 
an  object  in  its  hidden  region.  However,  what  is  the  dominant  factor  in  the  scat¬ 
tered  field?  Is  it  the  impedance  mismatch  at  r  =  b,  or  is  it  clue  to  a  scattered  field 
resulting  from  the  incident  field  being  transmitted  into  the  hidden  region  and  being 


reflected  by  the  hidden  object?  This  question  motivated  the  investigation  discussed 
in  the  following  paragraphs.  Initially,  this  work  attempted  to  minimize  the  energy 
transmitted  into  the  hidden  region.  It  was  theorized  less  energy  transmitted  into 
the  hidden  region  would  result  in  less  energy  that  can  be  scattered  by  the  cloaked 
object,  possibly  resulting  in  a  smaller  total  scattered  field.  The  constraints  on  the 
cylindrical  cloak’s  ideal  constitutive  parameters  defined  in  Equations  4.7  -  4.9  were 
used  to  define  parameter  sets  in  order  to  control  the  amount  of  energy  transmitted 
into  the  cloak’s  hidden  region.  The  cloaks’  effectiveness  was  then  analyzed  in  terms 
of  the  amount  of  energy  in  the  hidden  region  and  of  the  overall  scattering  width  of 
the  cloaking  structure. 

Simulations  using  Cornsol  Multiphysics  were  performed  on  the  original  simpli¬ 
fied  cylindrical  cloak  (Equation  4.5).  These  were  performed  with  a  PEC  cylinder 
with  radius  r  —  a  and  a  square  PEC  with  side  length  a  separately  in  the  cloak’s  hid¬ 
den  region.  These  objects  were  chosen  because  the  intent  was  to  show  objects  with 
different  scattering  properties  placed  in  the  hidden  region  impact  the  cloak’s  overall 
scattered  field.  Since  the  simulation  wavelength  is  A  =  a,  the  objects  placed  in  the 
hidden  region  are  on  the  order  of  one  wavelength  in  size.  The  chosen  objects  have 
significantly  different  shapes  and  areas;  thus,  the  overall  scattered  fields  should  differ 
due  to  field  penetration  into  the  hidden  region.  Additionally,  a  simulation  was  done 
with  no  objects  in  the  hidden  region.  All  simulation  results  are  shown  in  Figure  4.1. 
Note  the  fields  in  the  hidden  region  for  the  empty  cloak,  an  expected  result  based 
on  the  work  done  in  [101].  Based  on  these  images,  it  is  difficult  to  fully  comprehend 
the  size  and  pattern  of  the  scattered  field  for  each  geometry.  Therefore,  the  Cornsol 
simulation  results  were  transformed  to  a  far  zone  two-dimensional  scattering  width. 
The  scattering  width  for  each  geometry  is  plotted  as  a  function  of  6  in  Figure  4.2  (6  is 
the  bistatic  angle  with  9  =  0°  being  the  forward  scatter  direction).  Every  scattering 
width  plot  in  this  section  is  normalized  by  the  maximum  scattering  width  value  for 
an  uncloaked  PEC  cylinder  of  radius  a. 
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Figure  4.1:  Scattered  electric  field  magnitudes  for  a  simplified  cylindrical  cloak 

that  has  (a)  nothing,  (b)  a  PEC  cylinder  with  radius  a,  and  (c)  a  square  PEC  of  side 
length  a  in  the  hidden  region. 


Figure  4.2:  Normalized  scattering  width  for  a  simplified  cloak  with  various  hidden 
objects. 

Note  the  cloaked  PEC  cylinder  (red  line  in  Figure  4.2)  does  have  a  smaller 
scattering  width  than  an  uncloaked  PEC  cylinder  (blue  line).  Also  note  the  variation 
in  the  scattering  widths  between  the  cloaked  cylinder  and  square  (black  line).  This 
variation  is  due  to  different  objects  being  placed  in  the  hidden  region.  To  better 
see  how  the  scattering  width  is  changed  when  different  objects  are  inserted  into  the 
hidden  region,  the  difference  between  the  scattering  widths  of  a  cloaked  PEC  cylinder 
and  cloaked  PEC  square  were  plotted.  The  results  are  shown  in  Figure  4.3.  Obviously, 
changing  the  object  in  the  hidden  region  has  an  impact  on  the  overall  scattered  held 
for  this  set  of  constitutive  parameters.  The  question  is,  will  minimizing  the  held 
transmitted  into  the  hidden  region  reduce  the  overall  scattered  held  variations  caused 
by  the  different  hidden  objects? 
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Figure  4.3:  Scattering  width  difference  for  a  simplified  cloak  with  a  PEC  cylinder 
and  a  PEC  square  in  the  hidden  region. 


As  discussed  in  Section  2.2,  there  have  been  suggested  improvements  to  the 
original  simplified  constitutive  parameters  [14,101].  The  improved  set  of  constitutive 
parameters  for  TM~  incident  fields  put  forth  in  [102]  are  shown  in  Equation  4.12. 


flr  — 


(4.12) 


The  improved  parameter  set  was  developed  with  the  goal  of  reducing  the  overall 
scattering  width  of  the  cloaking  structure  by  matching  the  cloak’s  impedance  to  free 
space  at  r  =  b  while  still  satisfying  the  requirements  shown  in  Equations  4.7  and  4.8. 
As  with  the  initial  simplified  parameters,  the  improved  set  does  not  satisfy  the  third 
constraint  equation  shown  in  Equation  4.9.  Cornsol  simulations  of  a  cylindrical  cloak 
with  the  constitutive  parameters  shown  in  Equation  4.12  were  performed.  The  cloak 
was  simulated  with  no  object  in  the  hidden  region,  and  with  a  PEC  cylinder  of  radius 
a  and  a  square  PEC  with  side  length  a  separately  in  the  hidden  region.  The  results 
are  shown  in  Figure  4.4.  It  is  obvious  the  scattered  field  magnitudes  in  the  region 
r  >  b  are  larger  for  the  two  cloaks  with  objects  in  the  hidden  region.  To  better  show 
this,  the  scattered  field  results  were  transformed  to  the  far  zone.  These  results  are 
shown  in  Figure  4.5.  The  goal  of  reducing  the  overall  scattering  width  from  a  cloak 
with  simplified  parameters  was  achieved  using  the  constitutive  parameters  shown  in 
Equation  4.12.  This  is  due  to  the  matched  impedance  at  r  =  b.  However,  notice 
how  the  scattered  fields  have  a  greater  change  in  magnitude  when  the  objects  in  the 
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Figure  4.4:  Scattered  electric  field  magnitude  for  an  improved  cylindrical  cloak  that 
has  (a)  nothing,  (b)  PEC  cylinder  of  radius  a,  and  (c)  square  PEC  of  side  length  a 
in  the  hidden  region. 


Figure  4.5:  Normalized  scattering  widths  from  an  improved  simplified  cloak. 


hidden  region  are  changed.  As  before,  this  can  seen  by  comparing  the  difference 
in  the  scattering  widths  for  the  cloaked  PEC  cylinder  (red  line  in  Figure  4.5)  and 
the  cloaked  square  PEC  (black  line).  This  is  shown  in  Figure  4.6.  The  average 
difference  in  scattering  width  for  a  simplified  cylindrical  cloak  with  the  parameters 
shown  in  Equation  4.5  having  a  PEC  cylinder  and  square  in  the  hidden  region  was 
1.77  dB.  The  average  difference  in  scattering  width  for  the  improved  cylindrical  cloak 
with  the  parameters  shown  in  Equation  4.12  was  2.90  dB.  Therefore,  even  though 
overall  scattering  width  has  been  reduced  (due  to  the  matched  impedance  at  r  =  b), 
the  larger  variation  in  the  scattered  fields  when  different  objects  were  placed  in  the 
hidden  region  suggests  more  energy  is  being  transmitted  into  the  hidden  region  of  the 
cloak  with  the  improved  constitutive  parameters  than  the  hidden  region  for  the  cloak 
with  the  original  simplified  parameters.  This  implies  a  matched  impedance  at  r  =  b  is 
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Figure  4.6:  Scattering  width  difference  for  an  improved  simplified  cloak  with  a  PEC 
cylinder  and  a  square  PEC  in  the  hidden  region. 

more  important  for  signature  width  reduction  than  minimizing  the  held  transmitted 
into  the  hidden  region.  This  is  further  tested  in  the  next  section.  It  will  be  shown 
for  these  two  cloaks  that  the  value  of  fig  determines  the  size  of  the  held  transmitted 
into  the  hidden  region. 

4-2  Reducing  field  transmission  into  the  hidden  region 

As  previously  mentioned,  the  constitutive  parameters  of  an  ideal  cloak  for  TMZ 
incident  waves  must  satisfy  the  constraints  shown  in  Equations  4.7  -  4.9.  The  sim¬ 
plified  cloaks  in  the  literature  focus  on  satisfying  Equations  4.7  and  4.8.  Equation 
4.9  has  never  before  been  discussed,  likely  due  to  the  assumptions  used  when  the 
initial  set  of  simplified  parameters  was  put  forth.  In  what  follows,  the  importance 
of  Equation  4.9  is  analyzed  in  terms  of  overall  scattering  width  and  of  how  well  the 
hidden  region  is  shielded  from  incident  energy. 

If  one  hrst  assumes  a  cloak’s  constitutive  parameters  satisfy  Equations  4.7  and 
4.8,  Equation  4.9  can  be  written  in  a  more  compact  form.  This  is  shown  in  Equation 
4.13. 

he  T  ho — 7 - r  —  0  (4.13) 

r  (r  —  a) 

Note  that  a  <  r  <  b.  The  analysis  is  initially  confined  to  cloaks  with  a  constant 
value  for  fig  (i.e.  fig  —  0).  This  means  the  smaller  fig ,  the  less  error  there  will  be  for 
any  value  of  r  when  trying  to  satisfy  Equation  4.13.  The  left-hand  side  of  Equation 
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4.13  can  be  calculated  using  the  values  for  fig  for  the  initial  simplified  parameter 
set  (jio  =  1)  and  the  improved  parameter  set  (fig  =  ^).  These  are  plotted  as  a 
function  of  r  in  Figure  4.7.  For  this  plot,  a  =  1  and  6  =  2.  Two  additional  plots 
are  shown  in  this  graph,  and  these  will  be  discussed  later.  Since  the  ideal  value  for 


Figure  4.7:  The  calculated  values  for  the  left-hand  side  of  Equation  4.13  using  non¬ 
ideal  values  for  fig.  Values  for  fig  are  1  (red  line),  ^  (green  line),  (cyan  line), 
and  (^)3  (black  line). 

Equation  4.13  is  zero  for  all  values  of  r.  A  larger  calculated  value  for  the  left-hand 
side  of  Equation  4.13  using  a  non-ideal  value  for  fig  means  a  larger  deviation  in  the 
material  parameter  from  that  of  the  ideal  cloak.  Based  on  this  graph,  one  would 
expect  to  find  the  held  transmitted  into  the  hidden  region  for  the  cloak  with  the 
initial  simplified  parameter  set  ( fig  =  1)  to  be  less  than  the  hidden  region  held  for 
the  improved  parameter  set  ( fig  =  tzt)-  This  is  due  to  the  fact  the  value  for  fig  for 
the  initial  parameter  set  results  in  a  smaller  error  in  Equation  4.13  than  the  value  for 
fig  in  the  improved  parameter  set.  Hence,  changing  objects  in  the  hidden  region  for 
the  cloak  with  the  original  simplified  parameters  (Figure  4.1)  would  have  less  impact 
on  the  overall  scattered  held  than  changing  objects  in  the  hidden  region  of  a  cloak 
with  the  improved  parameters  (Figure  4.4).  This  is  precisely  what  we  have  shown  in 
Figures  4.3  and  4.6. 

Other  parameter  sets  can  satisfy  Equations  4.7  and  4.8  but  also  further  reduce 
the  deviation  from  Equation  4.13,  resulting  in  a  reduction  in  the  held  transmitted  into 
the  hidden  region.  The  parameter  sets  shown  in  Equations  4.14  and  4.15  meet  these 
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conditions.  Their  deviations  from  the  ideal  value  of  Equation  4.13  were  calculated. 
The  results  are  plotted  in  Figure  4.7. 
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Obviously  no  attempt  was  made  to  match  impedances  at  the  r  =  b  interface  for 
the  parameter  sets  shown  in  Equations  4.14  and  4.15,  as  the  goal  was  to  show  a 
reduction  in  the  field  transmitted  into  the  hidden  region.  The  cloaks  were  simulated 
with  material  parameters  shown  in  Equations  4.14  and  4.15  with  a  PEC  cylinder 
of  radius  a  and  with  a  square  PEC  with  side  length  a  in  the  hidden  region.  The 
difference  in  scattering  widths  are  plotted  for  each  of  these  cloaks.  The  results  are 
shown  in  Figures  4.8  and  4.9.  As  expected,  the  average  difference  in  scattering  width 


Figure  4.8:  Scattering  width  difference  for  a  cloak  with  parameters  shown  in  Equa¬ 
tion  4.14  with  a  PEC  cylinder  and  square  in  the  hidden  region. 


is  decreased  when  He  =  ^  and  He  =  (^)3  respectively,  leading  to  the  conclusion 
that  the  field  transmitted  into  the  hidden  region  is  being  reduced  as  He  is  decreased. 

To  further  show  how  Equation  4.13  determines  the  amount  of  energy  trans¬ 
mitted  into  a  simplified  cylindrical  cloak’s  hidden  region,  cloaks  with  the  material 
parameters  shown  in  Equations  4.5,  4.12,  4.14,  and  4.15  were  simulated.  No  objects 
were  placed  in  their  hidden  regions,  and  the  total  electric  field  magnitudes  in  the  hid¬ 
den  regions  for  each  cloak  are  plotted  in  Figure  4.10.  The  images  in  Figure  4.10  clearly 
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Figure  4.9:  Scattering  width  difference  for  a  cloak  with  parameters  shown  in  Equa¬ 
tion  4.15  with  a  PEC  cylinder  and  square  in  the  hidden  region. 

Table  4.1:  Hidden  Region  Total  Energy  and  Impedance  for  Different  Cloaks 


He 

Total  Energy 

Z\r=b 

b 

b—a 

2.76  pJ 

1 

l 

2.24  pJ 

0.5 

b—a 

b 

1.35  pJ 

0.25 

TXf 

\  b—a  ) 

0.42  pJ 

0.0625 

show  as  the  cloak  takes  on  values  of  fig  which  make  the  left-hand  side  of  Equation  4.13 
closer  to  zero  for  all  values  of  r,  there  is  less  field  transmitted  into  the  hidden  region. 
As  additional  proof,  the  energy  density  in  the  hidden  region  can  be  integrated  to 
determine  the  regions’  total  energies.  These  results  are  shown  in  Table  4.1.  Note  the 
cloak  with  material  parameters  shown  in  Equation  4.15  has  the  smallest  total  energy 
in  the  hidden  region.  Therefore  it  is  the  best  of  the  four  cloaks  considered  at  shield¬ 
ing  the  hidden  region.  However,  there  is  a  price  to  pay  for  this  improved  shielding 
performance.  Table  4.1  also  shows  the  impedance  for  each  cloak  at  r  =  b.  The  cloak 
with  the  best  shielding  of  the  hidden  region  also  has  the  worst  impedance  mismatch 
at  the  cloak  outer  boundary.  As  demonstrated  earlier,  an  impedance  mismatch  at  the 
boundary  results  in  the  cloaking  body  having  a  significant  scattered  held.  Therefore, 
to  compare  performance  in  terms  of  scattering  width  reduction,  scattering  widths 
were  plotted  for  all  cloaks  with  a  PEC  cylinder  of  radius  a  in  the  hidden  region. 
This  is  shown  in  Figure  4.11.  The  blue  line  is  the  normalized  scattering  width  for 
an  uncloaked  PEC,  the  red  line  is  the  normalized  scattering  width  for  cloak  with 
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Figure  4.10:  Electric  field  magnitude  in  the  hidden  region  for  cloaks  with  material 
parameters  defined  by  (a)  Equation  4.12,  (b)  Equation  4.5,  (c)  Equation  4.14,  and 
(d)  Equation  4.15. 


the  simplified  parameter  set  (Equation  4.5),  the  green  line  is  the  scattering  width  for 
a  cloak  with  the  improved  parameter  set  (Equation  4.12),  and  the  cyan  and  black 
lines  are  the  scattering  widths  for  cloaks  with  parameter  sets  shown  in  Equations 
4.14  and  4.15  respectively.  By  altering  the  material  parameters  such  that  less  fields 
are  transmitted  into  the  hidden  region,  the  change  in  impedance  at  r  =  b  dramat¬ 
ically  increases  the  overall  scattering  width  of  the  cloaking  structure.  Obviously  if 
scattering  width  reduction  is  the  goal,  use  of  the  improved  simplified  parameter  set 
(Equation  4.12,  [102])  is  the  best  option  as  its  scattering  width  is  significantly  less 
than  an  uncloaked  PEC  cylinder.  The  two  new  cloaks  actually  have  larger  scattering 
widths  at  various  observation  angles  than  the  uncloaked  PEC,  making  them  a  bad 
choice  if  signature  reduction  is  desired.  However,  if  one  is  attempting  to  simply  shield 
an  object  from  incident  radiation,  then  the  use  of  Equation  4.9  becomes  important 
in  that  parameters  should  be  chosen  such  that  the  deviation  from  this  equation  is 
minimized.  One  may  say  shielding  can  easily  be  accomplished  using  a  PEC;  why  use 
a  modified  cloak  for  such  a  task?  A  PEC  does  act  as  a  suitable  barrier  for  a  large 


67 


CO  10 


o 


z  -50 


0 


20 


40 


60 


80 


100 


120 


140 


160 


180 


6 


Figure  4.11:  Scattering  widths  for  cloaks  with  a  PEC  cylinder  of  radius  a  in  the 

hidden  region. 

bandwidth  of  electromagnetic  radiation.  However,  at  extremely  low  frequencies,  skin 
depths  must  be  considered.  It  might  be  less  costly,  in  terms  of  weight  or  size,  to  use 
a  designed  cloak  for  such  a  shielding  application. 

It  has  been  shown  that,  in  terms  of  overall  signature  reduction,  the  cloak  put 
forth  in  [102]  is  the  best  option,  and  that  this  particular  cloak  satisfies  Equations 
4.7  and  4.8  and  has  a  matched  impedance  at  r  =  b.  However,  note  the  significant 
variation  in  the  scattered  held  when  different  objects  are  placed  in  this  cloak’s  hidden 
region.  Is  it  possible  to  reduce  this  variation  in  the  scattered  held  with  different 
objects  in  the  hidden  region  while  maintaining  the  overall  reduction  in  scattering 
width? 

There  are  four  parameters  that  must  be  met  for  a  cylindrical  cloak  to  be  perfect: 
Equations  4.7,  4.8,  and  4.9  must  be  satished,  and  the  cloak  must  have  a  matched 
impedance  at  r  =  b.  To  this  point,  simplified  cloaks  that  satisfy  Equations  4.7  and 
4.8,  and  that  either  do  or  do  not  have  a  matched  impedance  at  r  =  b  have  been 
analyzed.  The  results  showed  the  cloak  with  the  matched  impedance  results  in  the 
best  improvement  in  scattering  width  even  though  this  cloak  has  the  largest  held 
transmitted  into  its  hidden  region. 

Now,  consider  a  cloak  with  material  parameters  shown  in  Equation  4.16. 


r  —  a 


r 


(4.16) 


68 


Like  the  improved  parameter  set  put  forth  in  [102],  these  parameters  satisfy  three  of 
the  four  requirements.  The  difference  is  these  parameters  satisfy  Equations  4.7  and 
4.9  while  having  a  matched  impedance  at  r  =  b.  Equation  4.8  is  not  satisfied. 

As  done  previously,  a  cloak  having  the  constitutive  parameters  shown  in  Equa¬ 
tion  4.16  was  simulated  with  a  PEC  cylinder  of  radius  a  and  a  square  PEC  of  side 
length  a  separately  in  the  hidden  region  were  simulated.  The  scattered  held  results 
are  shown  in  Figure  4.12.  Note  how  the  scattered  fields  for  all  three  images  shown  in 
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Figure  4.12:  Scattered  electric  held  magnitude  for  a  cylindrical  cloak  with  param¬ 
eters  shown  in  Equation  4.16  that  has  (a)  nothing,  (b)  PEC  cylinder,  and  (c)  square 
PEC  in  the  hidden  region. 


Figure  4.12  appear  very  similar.  This  suggests  different  objects  in  the  hidden  region 
have  little  effect  on  the  scattered  held.  This  can  seen  more  clearly  by  transforming 
the  scattered  helds  to  the  far  zone.  This  is  shown  in  Figure  4.13.  The  blue  line  is 


Figure  4.13:  Normalized  scattering  width  from  cloak  defined  by  Equation  4.16. 


the  scattering  width  for  an  uncloaked  PEC,  the  red  and  green  lines  are  the  scattering 
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widths  for  the  improved  cloak  with  a  PEC  cylinder  and  square  PEC  in  the  hidden  re¬ 
gion.  The  cyan  and  black  lines  are  the  same  but  for  a  cloak  with  material  parameters 
shown  in  Equation  4.16.  The  graphs  in  Figure  4.13  clearly  show  there  is  no  difference 
in  the  scattering  width  when  different  objects  are  placed  in  the  hidden  region  of  the 
cloak  with  parameters  put  forth  in  Equation  4.16.  Note  the  cyan  and  black  lines  lie 
virtually  on  top  of  the  other.  In  fact,  the  average  difference  in  the  scattering  widths 
is  0.0671  dB.  This  is  because  the  total  held  in  the  hidden  region  is  negligible.  The 
held  in  the  hidden  region  is  negligible  because  the  impedance  at  r  =  a  — >  oo,  which 
means  no  energy  will  be  transmitted. 

Figure  4.13  shows  scattering  width  results  for  a  cloak  using  the  parameters 
shown  in  Equation  4.12.  This  is  done  to  compare  the  performance  of  the  cloaks 
in  terms  of  scattering  width.  Obviously,  the  red  and  green  lines  are  more  desirable 
results  because  of  the  smaller  scattering  width  values.  However,  this  cloak’s  scattering 
widths  vary  more  as  a  function  of  different  objects  in  its  hidden  region.  If  various 
objects  are  going  to  be  hidden  and  observation  angles  are  in  the  specular  region,  the 
top  cloak  may  be  a  better  option.  Of  course,  the  cloak  with  parameters  put  forth 
in  Equation  4.16  has  two  radially  varying  parameters,  meaning  it  is  currently  more 
difficult  to  manufacture. 

Thus  far,  cylindrical  cloaks  that  satisfy  the  ideal  values  for  £z/ie  and  £z/j,r  have 
been  analyzed.  It  has  been  shown  deviations  from  a  third  constraint  equation  shown 
in  Equation  4.9  resulted  in  larger  fields  being  transmitted  into  a  cylindrical  cloak’s 
hidden  region.  As  the  cloak’s  constitutive  parameters  were  changed  such  that  this 
new  constraint  was  better  satisfied,  the  amount  of  energy  transmitted  into  the  hidden 
region  was  shown  to  be  reduced.  However,  the  resulting  impedance  mismatch  at  r  =  b 
due  to  changing  the  constitutive  parameters  resulted  in  a  significant  scattered  held. 
Thus,  despite  reducing  energy  transmitted  into  the  hidden  region,  which  resulted  in 
a  reduction  in  the  scattered  held  by  the  cloaked  object,  the  cloak  itself  was  creating 
a  large  scattered  held.  Hence,  in  terms  of  overall  scattering  width,  having  a  matched 
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impedance  at  r  —  b  seems  to  be  more  important  than  reducing  the  transmitted  energy 
into  the  hidden  region. 

A  particular  cylindrical  cloak  that  satisfied  the  specific  values  for  ezye  and  y'g 
while  having  a  matched  impedance  at  r  =  b  has  also  been  analyzed.  For  observation 
angles  in  the  backscatter  region,  this  cloak  to  performed  quite  well  in  terms  of  scat¬ 
tering  width  as  the  scattered  field  was  independent  of  objects  placed  in  the  hidden 
region.  While  scattering  width  performance  was  not  on  the  same  level  as  the  cloak 
with  parameters  put  forth  in  [102],  the  independence  of  the  scattered  field  due  to 
different  objects  in  the  hidden  region  is  noteworthy. 

4-3  Improved  Cylindrical  Cloak  Parameters 

As  discussed  in  Section  2.2,  simplified  parameter  sets  are  necessary  in  order  to 
make  cloaks  more  manufacturable.  Due  to  a  derivation  error  in  one  of  the  first  papers 
on  simplified  cloaks,  most  simplified  parameter  sets  have  focused  on  the  material 
parameters  satisfying  Equations  4.7  and  4.8  only.  As  shown  in  Section  4.1,  a  third 
equation  (Equation  4.9)  is  required  in  order  to  fully  constrain  the  material  parameters. 
Recall  also  from  Section  4.1  that  if  Equation  4.7  is  assumed  true  for  a  given  material 
parameter  set,  Equation  4.9  can  be  rewritten  as 

ID  +  Rfl— 7 - v  =  0-  (4-17) 

r  (r  —  a) 

Equation  4.17  is  only  a  function  of  the  r-varying  parameter,  yo ,  and  its  first  derivative, 
fi'() .  The  solution  to  this  first  order  differential  equation  is  given  by 

Mr)  =  C^—.  (4.18) 

r  —  a 

Note  C  is  a  constant.  Initially  C  is  assumed  to  be  unity  as  this  corresponds  to  a  first 
order  transformation,  as  will  be  discussed  below.  Obviously,  Equation  4.18  matches 
the  form  of  ye  shown  in  Equation  1.21.  What  may  not  be  as  obvious  is  the  fact  that 
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Tabic  4.2:  Simplified  material  parameter  values  for  N  —  0, 1,  and  2. 
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the  ideal  value  of  ezne  is  inherently  satisfied  in  Equation  4.17.  Thus,  to  develop  a 
simplified  parameter  set,  one  should  first  obtain  an  approximate  value  for  fig  that  is  as 
close  as  possible  to  the  solution  to  Equation  4.17.  To  obtain  values  for  jig  that  better 
approximate  the  ideal  value,  the  solution  of  the  differential  equation  for  He  can  be 
expanded  using  a  Taylor  series  about  the  point  r  =  b  [26].  Approximations  are  then 
made  by  limiting  N  in  the  Taylor  series.  Thus,  the  expression  for  the  approximation 
of  /ig  is 

N  77,  /  7  \ 

=  (4.19) 

Til 

n= 0 

where  He(b)  is  the  nth  derivative  of  /jl g(r )  evaluated  at  r  =  b.  Once  the  expression 
for  the  approximation  to  jig  is  determined,  values  for  Hr  and  are  found  using 
the  remaining  constraints  defined  in  Equations  4.7  and  4.8.  Since  the  expansion  is 
performed  about  the  point  r  =  b ,  the  impedance  at  r  —  b  is  matched  to  that  of  free 
space.  The  calculated  material  parameters  using  the  first  term  (N  =  0),  the  first  two 
terms  (N  =  1),  and  the  first  three  terms  (N  =  2)  to  estimate  jig  are  shown  in  Table 
4.2.  For  N  =  0,  jig  =  r^.  Consequently,  the  material  values  are  those  of  the  improved 
cloak  discussed  in  Section  2.2  and  originally  put  forth  in  [102],  This  explains  why 
the  improved  cloak  was  shown  to  be  more  effective  than  the  quadratic  and  original 
simplified  cloaks.  Note  also  that  while  these  expressions  may  seem  complicated,  the 
material  parameters  themselves  are  quite  well  behaved.  This  is  shown  in  Figure  4.14. 
As  additional  terms  are  included  in  the  Taylor  series,  He  and  become  more  varying 
with  respect  to  r;  they  also  take  on  more  extreme  values  at  r  =  a.  Hence  these 
parameter  sets  can  be  used  based  on  the  ability  to  manufacture  metamaterials  with 
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N  =  1 


Figure  4.14:  Parameter  variation  as  a  function  of  r.  In  both  figures,  the  o-line  is 
he,  the  A-line  is  fir ,  and  the  D-line  is  ez.  In  the  top  figure  fie  is  approximated  using 
N  =  1;  in  the  bottom  figure,  N  =  10. 

various  properties.  If  larger  values  of  fig  are  attainable,  simply  let  N  increase  to  define 
a  new  parameter  set. 

The  analysis  up  to  this  point  has  assumed  a  first-order  transformation  as  shown 
in  Equation  1.7.  It  can  be  expanded  to  include  all  transformation  orders.  As  an 
example,  an  nth-order  coordinate  transformation  that  maps  the  region  r'  <  b  to  the 
region  a  <  r  <  b  has  the  form 


r  = 


1 

bn~i 


¥]r'n  +  a- 


(4.20) 


The  material  parameters  for  a  TAP  electromagnetic  field  can  be  found  using  the 
method  in  [72]  and  are  shown  below. 


r  —  a 

fir  =  n - 

r 


he 


1  r 
nr  —  a 


(r  —  a) n  1  b2 
rn  ( b-a Y 


(4.21) 


As  done  above,  the  ideal  material  parameters  shown  in  Equation  4.21  can  be  substi¬ 
tuted  into  the  general  wave  equation  shown  in  Equation  4.1  in  order  to  obtain  the 
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constraints  on  the  material  parameters. 


1  {b  —  a)n  n 2 

^zk’d  (j-  —  d'jri  2 


(4.22) 


1  (b  —  a) n  r2 

Sz^r  b2  (j.  —  (i) n 


i  i  iA 

£zl*e  r  £zy2d 


(b  —  a)n  n2 
u  (r  —  a)  n 


1  a 

r  r  (r  —  a) 


(4.23) 

(4.24) 


As  before,  Equation  4.24  can  be  simplified  if  Equations  4.22  is  assumed  to  be  true. 
The  result  is 

N  +  — r  =  0,  (4.25) 

r  (r  —  a ) 

which  is  identical  to  Equation  4.13.  Hence,  the  choice  of  C  in  Equation  4.18  is 
important  because  C  =  ^  where  n  is  the  transformation  order.  If  C  ^  1,  then  one 
cannot  use  Equations  4.7  and  4.8  when  solving  for  and  /ir .  Rather,  Equations  4.23 
and  4.24  must  be  used  to  get  the  proper  approximations  for  the  material  parameters. 


In  the  next  section,  the  performance  of  cloaks  where  no  is  approximated  using 
different  values  of  N  in  the  Taylor  series  are  compared. 


4-4  Analysis  of  Cloak  Performance 

The  Cornsol  Multiphysics  software  package  was  used  to  perform  all  simulations. 
Cylindrical  cloaks  whose  material  parameters  were  found  using  the  process  described 
in  Section  4.3  using  a  linear  transformation  were  simulated.  The  cloaks’  inner  bound¬ 
aries  were  lined  with  a  PEC  and  an  incident  wavelength  (A)  of  one  meter  traveling  in 
the  positive  x  direction  was  used.  The  cloak  parameters,  a  and  b ,  were  dehned  such 
that  a  —  A  and  b  =  2A.  The  impact  on  the  results  when  a  and  b  were  varied  is  dis¬ 
cussed  later.  Simulation  results  for  the  improved  cloak  and  that  of  a  cloak  with  a  10, 
50,  and  100-term  Taylor  series  approximations  for  /ig  are  shown  in  Figure  4.15.  Note 
all  cloaks  in  Figure  4.15  show  good  cloaking  performance.  However,  it  is  difficult  to 
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1.246 


Figure  4.15:  Real  part  of  the  5-component  of  the  total  electric  field  for  cloaks  with 
material  parameters  defined  as  (a)  the  improved  cloak,  while  (b),  (c),  and  (d)  are 
cloaks  where  fig  is  approximated  using  a  10,  50,  and  100-term  Taylor  series. 

determine  which  cloak  has  the  best  performance.  Therefore,  only  the  scattered  field 
outside  of  the  cloaks  for  each  configuration  were  plotted.  These  results  are  shown 
in  Figure  4.16.  Figure  4.17  shows  the  scattered  fields  from  cloaks  with  the  same 
configurations  except  that  a  second  order  transformation  was  used  when  developing 
the  ideal  parameter  values.  Based  on  the  analysis  in  Section  4.3,  the  improved  cloak 
should  have  the  largest  scattered  field  while  the  cloak  where  fig  is  approximated  using 
a  100-term  Taylor  series  should  have  the  smallest  scattered  field  for  all  transformation 
orders.  This  is  clearly  evident  in  both  figures,  particularly  when  one  compares  the 
forward  scattered  fields. 

Based  on  these  results,  one  can  conclude  the  process  discussed  in  Section  4.3 
results  in  material  parameter  sets  that  result  in  better  cloak  performance  as  the  num¬ 
ber  of  terms  in  the  Taylor  series  is  increased  for  all  transformation  orders.  However, 
in  order  to  better  classify  cloak  performance,  the  near-field  results  were  transformed 
to  the  far-held  in  order  to  determine  each  cloak’s  two-dimensional  scattering  width. 
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Figure  4.16:  Real  part  of  the  ^-component  of  the  scattered  electric  held  for  cloaks 
with  material  parameters  defined  as  (a)  the  improved  cloak,  while  (b),  (c),  and  (d) 
are  cloaks  where  /x#  is  approximated  using  a  10,  50,  and  100-term  Taylor  series. 

In  order  to  compare  the  cloak’s  scattering  width  reduction  capabilities,  a  met¬ 
ric  must  be  established.  First  the  best  reduction  in  scattering  width  numerically 
possible  was  determined.  This  occurs  when  an  ideal  cylindrical  cloak  is  used,  theo¬ 
retically  resulting  in  no  scattered  held.  However,  due  to  numerical  issues,  there  is  a 
residual  scattered  held.  Specifically,  the  finite  element  method  approximates  the  held 
behavior  with  piecewise  continuous  elements.  Therefore  this  discretization  limits  the 
accuracy  of  the  held  representation,  particularly  near  r  =  a  where  the  cloak  has  pa¬ 
rameter  values  equal  to  zero  or  infinity.  Approximations  must  be  made  because  these 
values  cannot  be  simulated  due  to  resulting  singularities  in  the  differential  equation. 
As  stated  in  Section  4.3,  analysis  has  shown  even  slight  deviations  from  the  ideal 
material  parameters  result  in  degraded  cloak  performance  [42,77].  To  illustrate  this 
degradation,  the  scattering  width  for  a  simulation  with  no  scattering  objects  in  the 
domain  was  plotted  and  compared  to  the  results  for  an  ideally  cloaked  PEC  cylinder. 
The  results  are  shown  in  Figure  4.18.  Note  all  plots  are  normalized  by  the  maxi¬ 
mum  scattered  held  magnitude  for  the  uncloaked  PEC  cylinder.  The  ideal  cloak  does 
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Figure  4.17:  Real  part  of  the  ^-component  of  the  scattered  electric  field  for  cloaks 
with  material  parameters  derived  using  a  second  order  transformation.  Images  (a), 
(b),  (c),  and  (d)  are  cloaks  where  hq  is  approximated  using  a  1,  10,  50,  and  100-term 
Taylor  series. 

significantly  reduce  the  scattering  width  of  the  uncloaked  PEC  cylinder.  However, 
one  would  expect  the  results  to  be  on  the  order  of  the  scattered  held  from  an  empty 
domain.  Clearly,  they  are  not  even  close  to  this  result  due  to  the  aforementioned 
approximations. 

The  simulation  can  be  improved  by  controlling  the  size  of  the  mesh  elements 
particularly  near  the  region  r  =  a.  A  new  simulation  domain  was  constructed  by 
creating  a  subdomain  around  r  =  a  of  extremely  small  elements.  This  increased  the 
resolution  with  regard  to  the  singular  parameter  values.  The  final  mesh  resulted  in 
approximately  1,500,000  degrees  of  freedom.  At  this  discretization,  the  ideal  results 
were  -30  dB  down  from  the  uncloaked  cylinder.  Further  improvements  required  sig¬ 
nificant  computation  times  and  would  not  further  aid  in  the  performance  analysis  as 
each  cloak  design  is  affected  equally.  Therefore  -30  dB  is  considered  to  be  the  ideal 
solution  for  comparison  purposes. 
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Figure  4.18:  Scattering  widths  for  an  uncloaked  PEC  cylinder  (o-line),  ideally 

cloaked  PEC  cylinder  (o-line),  and  an  empty  domain  (A-line). 

The  performance  of  the  improved  cloak  was  compared  to  cloaks  whose  mate¬ 
rial  parameters  were  found  using  the  process  described  in  Section  4.3.  Results  for 
the  improved  cloak  and  that  of  a  cloak  with  a  10,  50,  and  100-term  Taylor  series 
approximations  to  fig  are  shown  in  Figure  4.19.  Note  even  the  10-term  Taylor  series 
approximation  results  in  a  three  dB  scattering  width  improvement  in  the  region  near 
9  =  0°  compared  to  the  improved  cloak.  Also  note  the  scattering  width  gets  better 
as  more  terms  in  the  Taylor  series  are  used,  which  is  expected  since  the  material 
parameters  become  the  ideal  ones  as  N  — >  oo. 

In  the  backscattering  region  ( 9  =  180°),  improvement  towards  the  ideal  solution 
is  extremely  slow  due  to  the  large  slope  in  the  ideal  value  for  fig  as  r  — >  a.  The  Taylor 
series  requires  more  terms  in  order  to  accurately  represent  fig  in  this  area.  It  is  also 
interesting  to  note  these  new  material  parameters,  although  all  are  spatially  varying, 
are  very  manufacturable.  For  example,  at  r  —  a,  fig  —  11  and  e2  =  0.3636  for  the 
10-term  approximation.  As  the  number  of  terms  increases,  the  values  become  harder 
and  harder  to  obtain  from  a  manufacturing  perspective  (fig  =  51  and  ez  =  0.0784  for 
50-term  approximation,  fig  =  101  and  £2  =  0.0396  for  the  100-term  approximation). 
However,  as  advances  in  material  manufacturing  technology  continue,  these  values  are 
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Figure  4.19:  Scattering  widths  for  an  uncloaked  PEC  cylinder  (o-line),  and  cloaked 
PEC  cylinders  using  the  ideal  cloak  (A-line),  the  improved  cloak  (□-line),  and  10,  50, 
and  100  term  approximations  for  [1q  (o,  *,  and  V-lines). 

becoming  more  approachable.  Note  similar  results  are  obtained  for  far  held  patterns 
when  a  second  order  coordinate  transformation  is  used. 

As  stated  above,  all  simulations  were  performed  with  A  =  1  meter  with  the  cloak 
inner  and  outer  boundaries  located  at  a  =  A  and  b  =  2A  respectively,  ft  was  shown 
in  [102]  the  improved  cloak  is  relatively  immune  to  changes  in  cloak  thickness.  As 
shown  in  Section  4.3,  the  improved  cloak  parameters  can  be  derived  using  the  Taylor 
series  expansion  method  by  letting  IV  =  0.  The  material  parameter  sets  developed 
in  this  section  are  refinements  to  the  improved  cloak  and  are  derived  by  increasing 
N  in  the  Taylor  series.  Therefore,  these  new  material  parameter  sets  should  also  be 
relatively  immune  to  changes  in  b.  A  number  of  simulations  were  performed  with  b 
varying  from  1.5 A  to  4. 5 A.  In  all  simulations,  the  results  were  similar  to  those  shown 
in  Figure  4.19.  Scattering  width  results  for  the  geometry  when  b  =  4.5A  are  shown 
in  Figure  4.20.  Note  there  is  a  slight  degradation  in  performance  as  b  is  increased. 
When  b  =  2A,  the  scattering  width  at  6  =  0°  was  -16  dB.  When  b  is  increased  to 
4. 5 A,  the  scattering  width  increased  to  -14  dB.  Hence,  a  thicker  cloak  will  result  in 
some  degradation  of  performance  for  the  same  material  parameter  set.  However,  for 
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Figure  4.20:  Scattering  width  for  cloaks  where  b  =  4.5A.  Plots  show  an  uncloaked 
PEC  cylinder  (o-line),  an  ideally  cloaked  PEC  cylinder  (A-line),  and  cloaks  for  yo 
approximations  using  10,  50,  100,  and  150  terms  (□,  o,  *,  and  V-lines). 

all  values  of  b ,  as  the  number  of  terms  goes  to  infinity,  the  performance  approaches 
that  of  the  ideal  case,  as  expected. 

4  •  5  Summary 

In  this  section,  a  new  way  to  develop  simplified  material  parameter  sets  for 
cylindrical  cloaks  was  presented.  Specifically,  for  TM2  incident  waves,  the  approx¬ 
imation  of  yo  should  first  be  defined  using  a  Taylor  series  expansion  of  the  ideal 
parameter  as  defined  using  all  constraint  equations.  TE2  results  are  easily  obtained 
by  the  application  of  duality.  The  constitutive  parameters  /ir  and  ez  can  be  deter¬ 
mined  by  making  the  products  yo£z  and  yr£z  equal  to  the  same  products  using  the 
ideal  material  parameter  set.  The  performance  of  cloaks  developed  in  this  manner 
is  limited  only  by  the  number  of  terms  used  in  the  Taylor  series  expansion,  which  is 
dictated  by  existing  manufacturing  capabilities.  Scattering  width  improvement  was 
observed  for  all  angles  when  compared  to  previous  published  material  parameter  sets. 
Significant  improvement  was  noted  in  the  forward  scattering  region.  It  was  also  shown 
the  simplified  parameter  set  put  forth  in  [102]  is  a  simplification  of  this  method  in 
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which  the  Taylor  series  expansion  of  fig  is  limited  to  the  first  term.  These  parameter 
sets  also  have  relatively  consistent  performance  for  all  values  of  b.  Performance  for 
a  constant  number  of  terms  in  the  Taylor  series  does  slightly  degrade  as  b  increases, 
but  for  all  b,  ideal  cloaking  performance  is  approached  as  N  — >  oo. 
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V.  Computational  Improvement  Using  a  Green’s  Function 

As  discussed  in  Chapter  II,  most  of  the  simulations  involving  cloaks  have  been  done 
using  the  Cornsol  Multiphysics  software,  a  finite-element-based  software  package.  The 
finite  element  method  is  particularly  useful  to  simulate  cloaks  due  to  the  inhomogene¬ 
ity  of  the  constitutive  parameters.  However,  there  are  limitations  and  trade-offs  which 
must  be  considered  when  solving  a  problem  using  FEM.  FEM  partitions  a  computa¬ 
tional  domain  into  a  large  number  of  subdomains  of  simple  geometric  shapes.  The 
partitioning  grid  is  commonly  called  the  mesh.  For  two-dimensional  problems,  tri¬ 
angular  elements  are  commonly  used,  but  there  are  others.  The  solution  to  the  dif¬ 
ferential  equation  is  obtained  for  each  element.  Thus,  the  finer  the  mesh  used  when 
modeling  the  problem  geometry,  the  more  accurate  the  solution  [95].  Some  geometries 
require  greater  mesh  fidelity  due  to  a  rapid  variation  in  geometry  or  spatial  parame¬ 
ters.  Similarly,  simulations  with  large  computational  domains  require  more  elements 
due  to  the  problem  size.  Larger,  denser  meshes  result  in  an  increased  computational 
burden  which  can  result  in  long  solution  times.  In  this  chapter,  a  Green’s  function 
formulation  is  developed  to  determine  scattering  widths  from  a  cylindrically  cloaked 
PEC  cylinder.  Solution  time  is  significantly  reduced  when  solving  for  <j2 d  using  the 
Green’s  function  compared  to  FEM  methods.  The  tradeoff  accompanying  the  Green’s 
function  implementation  is  the  cloak  formulation  is  now  limited  to  a  cylindrical  ge¬ 
ometry  with  a  PEC  lining  the  inner  boundary.  Changing  the  cloak  geometry  would 
require  the  derivation  of  a  different  Green’s  function. 


5. 1  Solution  Geometry 

Consider  the  geometry  shown  in  Figure  5.1.  A  theoretical  solution  for  the  PEC 
cylinder’s  scattering  width  when  illuminated  by  an  incident  TEZ  plane  wave  exists 
and  has  the  closed  form  analytic  expression  shown  below  [13]. 


r  2A 
u  2D  =  hm  — 

r— >  oo  7T 


oo 


/  / 1 m  „r2V/, — 7  cos (mO) 

m— 0 


1  rr(2 )' 


(■ Ka ) 


(5.1) 


82 


Figure  5.1:  Geometry  for  an  FEM  simulation  of  scattering  from  a  PEC  cylinder. 

In  Equation  5.1,  kQ  is  the  free  space  wave  number,  6  is  the  observation  angle,  and  em 
is  1  for  m  =  0  and  2  otherwise.  Note  also  that  '  implies  differentiation  with  respect  to 
r.  Obviously  Equation  5.1  is  an  infinite  sum,  meaning  the  exact  theoretical  solution 
can  never  be  computed  using  a  computer.  However,  m  can  be  truncated  based  on  the 
specified  level  of  accuracy.  For  this  work,  the  summation  in  Equation  5.1  is  truncated 
to  M  using  the  following  criteria: 

x  =  max  \Fm\  ,  m  G  [0,  M\  , 

V  m  >  M,  \Fm\  <  O.Olx.  (5.2) 

where  Fm  =  J™(k°a'>  in  Equation  5.2.  The  validity  of  truncation  using  Equation  5.2 

{k0a) 

can  be  verified  by  comparing  the  calculated  scattering  widths  of  a  PEC  cylinder  of 
radius  a  =  X  for  rn  —  10,  which  is  the  determined  sum  limit  based  on  the  criteria  in 
Equation  5.2,  and  m  =  50.  The  metric  used  to  compare  the  similarity  between  the 
solutions  is  the  average  difference  in  <T2d-  Mathematically,  this  is  written  as 

1  M 

A  =  (<J2d(^p)  ~  °2 L(^p))  j  (5-3) 

P= 1 
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TEZ  Scattering  from  a  PEC  Cylinder 


Figure  5.2:  Scattering  widths  for  PEC  cylinder  calculated  using  analytic  and  FEM 
solutions. 

where  M  is  the  total  number  of  observation  angles,  and  the  <72 d  terms  are  the  scat¬ 
tering  widths  which  will  be  compared.  A  for  a^D  =  cr2D|m=io  and  cr =  cr2£>|m=50  is 
0.0063  m2.  The  scattering  width  for  m  =  10  is  plotted  in  Figure  5.2  (blue  line).  Note 
u 2d  is  on  the  order  of  1  -  20  m2.  Thus,  a  A  of  0.0063  m2  is  considered  negligible. 

The  analytic  solution  in  Equation  5.1  can  be  compared  with  the  FEM  solution 
obtained  using  the  Cornsol  software.  Note  the  computational  boundary  is  only  5A  x  5A. 
As  was  done  in  Chapter  IV,  the  near  held  results  obtained  using  Cornsol  can  be 
transformed  to  the  far  zone  to  obtain  the  two-dimensional  scattering  width.  Different 
meshes  were  used  in  the  Cornsol  simulations.  A  smaller  maximum  element  length 
(MEL)  corresponds  to  a  finer,  denser  mesh.  Simulation  results  were  obtained  using 
seven  different  meshes.  Results  are  plotted  in  Figure  5.2  with  additional  information 
on  problem  size  and  solution  speed  listed  in  Table  5.1. 

Obviously  there  is  good  visual  agreement  between  the  analytic  and  FEM  so¬ 
lutions  with  the  only  noticeable  differences  occurring  when  MEL  =  0.5A.  However, 
a  metric  was  needed  other  than  visual  alignment  to  define  FEM  solution  accuracy. 
Therefore,  based  on  the  results  in  Figure  5.2  and  Table  5.1,  A  on  the  order  of  0.1 
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Table  5.1:  Analytic  and  FEM  Solution  Comparisons 


MEL 

Unknowns 

Time 

Ao'2dL,=io 

0.5A 

2,056 

0.17  s 

0.171  m2 

0.25A 

3,832 

0.23  s 

0.125  m2 

0.1A 

23,146 

1.04  s 

0.039  m2 

0.075A 

40,654 

1.67  s 

0.039  m2 

0.05A 

92,168 

4.02  s 

0.043  m2 

0.025A 

366,024 

16.5  s 

0.043  m2 

0.01A 

2,292,872 

1,493  s 

0.043  m2 

m2  was  defined  as  the  threshold  for  good  solution  agreement.  For  the  geometry  in 
Figure  5.1,  an  MEL  of  0.1A  is  sufficient  to  obtain  good  solution  agreement.  Addi¬ 
tional  mesh  fidelity  does  not  result  in  better  solution  agreement,  and  it  also  requires 
significantly  longer  solution  times.  As  will  be  seen  shortly,  there  are  cloak  geometries 
where  an  MEL  =  0.1  A  does  not  result  in  sufficient  mesh  fidelity  due  to  the  thinness 
of  subdomains  within  the  computation  region. 


5.2  Green’s  Function  for  a  Layered  PEC  Cylinder 

Green’s  functions  are  routinely  used  in  electromagnetic  scattering  problems. 
However,  they  have  not  been  applied  to  solve  radiation  problems  involving  cloaks, 
likely  due  to  the  difficulty  in  their  derivation  due  to  the  anisotropic  nature  of  a 
cloak’s  material  parameters.  As  discussed  in  Section  2.4,  a  cylindrical  cloak  can 
be  approximated  by  using  concentric  layers  of  isotropic  material  with  homogeneous 
permittivity  and  permeability.  Therefore,  a  Green’s  function  for  a  magnetic  line 
source  in  the  far  zone  radiating  in  the  presence  of  a  PEC  cylinder  of  radius  r  —  a 
covered  by  n  layers  of  homogeneous  material  approximates  a  TE2  plane  wave  incident 
on  a  cloaked  cylinder.  This  geometry  is  shown  in  Figure  5.3.  Once  the  Green’s 
function  is  known,  the  total  field  can  be  found  by 

^  TOT  k2  <-  «-» 

H  =  ■  G,  (5.4) 

4:UJ  fJjQ 
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Magnetic  Line 
Source  in  Far  Zone 


Figure  5.3:  Problem  geometry  for  Green’s  function  derivation 


where  Im  is  the  strength  of  the  magnetic  line  source,  /  is  the  unit  dyad,  and  G  is  the 
Green’s  function. 

The  derivation  of  such  a  Green’s  function  is  not  a  new  concept.  Methods  similar 
to  what  is  proposed  here  have  been  used  to  study  near  and  far  held  solutions  for  a 
PEC  cylinder  covered  by  isotropic  lossless  materials  [10,103].  However,  these  analyses 
focused  on  the  scattering  properties  of  PEC  cylinders  layered  with  double  negative 
materials. 

A  Green’s  function  for  the  geometry  shown  in  Figure  5.3  was  derived  using  the 
method  described  in  [23]  with  the  details  shown  in  Appendix  B.  A  Green’s  function 
for  a  magnetic  line  source  radiating  in  the  presence  of  a  PEC  cylinder  with  n  layers 
of  homogeneous,  isotropic  material  has  the  form: 


G  =  ^  em  cos  [m{9  -  9')} 


m= 0 


Bn+1 

MKr<)  +  0r<) 


flt2)(*vr>),  (5-5) 
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(5.6) 

(5.7) 


Note  that  r<,  r>  are  the  lesser  and  greater  of  r  and  r'  respectively,  and  9 ,  9’  are  the 
observation  and  source  angular  locations.  The  remaining  unknowns  in  Equation  5.5 


are  the  B ^+1  coefficients.  These  coefficients  are  determined  based  on  the  junction 
conditions  at  the  radial  boundaries  which  force  continuity  of  tangential  magnetic  and 
electric  fields.  Due  to  the  PEC  boundary  at  r  =  a,  the  Bxm  value  is  known,  which 
allows  the  remaining  values  to  be  found  by  solving  a  system  of  equations  of  the  form 
Ax  =  B,  where  A  is  a  2 n  x  2 n  matrix,  and  n  is  the  number  of  layers  surrounding  the 
PEC  (see  Appendix  B  for  details).  Note  Equation  5.5  is  valid  when  observing  the 
field  where  r  >  rn  i.e.  in  the  free  space  region. 

Equation  5.5  contains  components  for  the  incident  field  and  the  scattered  field. 
The  incident  field  is  represented  by  the  Jm{k0r< )  component  while  the  scattered  field 
is  represented  by  the  i/m')(A;0r<)  terms.  Thus,  the  Green’s  function  can  be  rewritten  as 
two  separate  functions.  This  is  done  because  the  goal  is  to  compute  d  and  compare 
the  result  obtained  using  a  Green’s  function  formulation  to  the  result  obtained  using 
an  FEM-based  method  (Cornsol).  Also,  without  loss  of  generality,  the  magnetic  line 
source  is  assumed  to  be  at  9'  =  180°  in  the  far  zone,  which  results  in  a  plane  wave 
traveling  in  the  x  direction.  Hence,  the  Green’s  function  for  the  incident  field  can 
first  be  written  as 


.  VXJ 

G*  =  Ee-(-l)mc°s  (m0)Jm(kor)H%\kor').  (5.8) 

m= 0 


The  large  argument  approximation  for  the  Hankel  function  [7]  can  then  be  used  to 
simplify  the  expression.  This  is  done  because  the  line  source  is  assumed  to  be  in  the 
far  zone  in  order  for  plane  wave  incidence. 


H£>(kor') 
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VT  kn 
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(5.9) 


The  following  is  defined  for  simplicity  of  writing. 


H0  = 


2  e 


-jkar' 


nkn 


-e3  4 


(5.10) 
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Hence,  the  Green’s  function  for  the  incident  held  at  the  surface  of  the  layered  PEC 
cylinder  (r  =  b)  is 


G'  =  - ]-H0  cos  (mO)Jm(kc b). 


(5.11) 


m= 0 


Similarly,  the  Green’s  function  for  the  scattered  held  can  be  represented  as 


Bl+1 


Gs  =  -jH0J2U-j)mcos(me)-^H%\k0r).  (5.12) 


m= 0 


The  large  argument  approximation  for  the  Hankel  function  can  again  be  applied  since 
the  ultimate  goal  is  to  calculate  a2D. 


i  I  2  p~ik°r  ,  Bn+l 

=  -  W  s;  -jre  ‘  E  cos  (m9)  m 

'  v  m= 0 


/I  71+1 


(5.13) 


In  general,  scattering  width  can  be  found  by 

\jjs\ ^ 

<y 2d  =  hm  27T r- — -o-  (5-14) 

r— >oo  |  | z 

Therefore,  the  scattering  width  for  a  layered  PEC  cylinder  with  an  incident  TE2  plane 
wave  traveling  in  the  positive  x  direction  is 
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(5.15) 


The  summation  in  Equation  5.15  is  truncated  to  m  =  M  using  the  criteria  put  forth 
in  Equation  5.2.  The  closed  form  analytic  solution  for  (t2d  can  be  used  to  solve  for 
the  scattering  width  from  a  layered  PEC  cylinder. 


5.3  FEM  and  Green’s  Function  Comparison 


Since  it  is  not  yet  possible  to  manufacture  cloaks  with  the  required  spatially 
varying  anisotropic  parameters,  concentric  rings  of  anisotropic  material  are  used  to 
approximate  the  spatial  variation  [79].  Such  a  realization  can  be  simulated  in  Comsol, 
but  the  thinness  of  the  layers  necessitates  more  elements  resulting  in  an  increased  com¬ 
putational  burden,  as  will  be  shown  below.  Two  thin  layers  of  homogeneous  isotropic 
material  can  be  used  to  approximate  the  anisotropic  concentric  rings.  This  geometry 
approximation  can  be  quickly  solved  using  the  Green’s  function  implementation. 

Consider  a  two-dimensional  ideal  cylindrical  cloak  with  material  parameters, 
repeated  below  for  convenience. 

r  -  a  r  r  -  a  f  b  \2 

£r  = - ,  ee  = - ,  Hz  = -  - -  (5-16) 

r  r  —  a  r  \b  —  a  J 

In  order  to  manufacture  this  cloak,  concentric  rings  of  anisotropic  material  would  be 
required  to  approximate  the  spatial  variation,  much  like  what  was  done  in  [79].  Based 
on  the  results  in  [40],  the  anisotropic  concentric  rings  can  be  approximated  by  using 
thin  layers  of  homogeneous  isotropic  material.  Recall  from  Section  2.4,  the  required 
material  parameters  for  the  layers  approximating  the  anisotropic  layer  are 


£e 


e 1  +  y£2 

1  +  Tj 


(5.17) 


±  =  *  (i  +  T), 

£r  1  +  ?7  \£i  £2/ 


(5.18) 


As  an  example,  consider  an  ideal  cloak  realized  using  ten  layers  of  concentric 
anisotropic  rings.  Twenty  layers  of  homogeneous,  isotropic  rings  can  be  used  to  sim¬ 
ulate  the  ten-layer  anisotropic  cloak.  The  values  of  £r ,  £g ,  and  pz  for  each  anisotropic 
ring  are  determined  by  evaluating  Equation  5.16  at  each  of  the  ten  layers  (r  =  rn ). 
The  permittivity  values  for  the  isotropic  layers  are  then  determined  by  substituting 


FEM  and  GF  c>2D  Comparison,  20  Layers 


Figure  5.4:  Green’s  function  and  FEM  results  comparison  for  PEC  with  20  layers 

the  values  for  er  and  £g  into  Equations  5.17  and  5.18.  The  value  of  /iz  is  determined 
by  evaluating  /iz  given  in  Equation  5.16  at  r  =  rn. 

A  20-layer  isotropic  approximation  of  the  ideal  cloak  surrounding  a  PEC  cylin¬ 
der  was  simulated  using  Comsol.  Scattering  width  values  were  determined  based  on 
Comsol  results.  Analytic  results  were  found  using  Equation  5.15.  All  results  are 
shown  in  Figure  5.4.  Additionally,  results  from  the  two  methods  for  a  simulation 
using  40  layers  of  isotropic  material  to  simulate  a  20-layer  cloak  of  anisotropic  con¬ 
centric  rings  approximating  an  ideal  cloak  were  computed.  These  are  shown  in  Figure 
5.5. 

Note  the  similarities  between  the  Comsol  and  Green’s  function  solutions.  For 
the  20-layer  results,  A  was  0.14  m2,  which  shows  good  agreement  between  the  two 
solutions.  However,  there  is  a  noticeable  difference  in  computation  time.  The  Green’s 
function  solution  took  2.28  s.  The  Comsol  solution  was  obtained  by  first  creating  a 
non-uniform  mesh  over  the  computational  domain.  This  was  necessary  because  the 
spacing  between  layers  was  only  0.05A  and  0.025A  for  the  20  and  40-layer  simulations 
respectively,  and  a  uniform  mesh  with  MEL  <  0.01A  was  not  possible  due  to  memory 
limitations  of  the  Dell  690  precision  work  station  with  eight  gigabytes  of  RAM.  The 
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FEM  and  GF  c>2D  Comparison,  40  Layers 


Figure  5.5:  Green’s  Function  and  FEM  Results  Comparison  for  PEC  with  40  Layers 

MEL  in  the  concentric  layers  was  limited  to  0.007A,  while  for  the  rest  of  the  com¬ 
putational  domain,  MEL  =  0.05A.  For  the  20-layer  FEM  simulation,  the  mesh  had 
911,004  elements  and  a  solution  time  of  125  s. 

The  40-layer  simulation  resulted  in  a  A  of  0.04  m2.  The  Green’s  function 
solution  time  was  2.82  s.  The  Comsol  solution  used  a  similar  non-uniform  mesh 
with  881,892  elements  and  had  a  solution  time  of  124  s.  Obviously,  the  Green’s 
function  method  is  faster,  particularly  if  a  number  of  simulations  are  to  be  performed 
to  conduct  an  optimization  or  an  error  analysis  based  on  parameter  or  thickness 
variations  in  the  layers.  Additionally,  if  more  layers  are  to  be  used,  an  FEM  solution 
will  require  finer  meshing  within  the  layers,  increasing  the  number  of  unknowns  which 
will  increase  solution  time. 

Another  benefit  of  using  the  Green’s  function  to  calculate  scattering  widths  is 
when  large  cloak  geometries  are  simulated.  Up  until  this  point,  all  previous  simula¬ 
tions  in  this  section  have  used  the  cloak  parameters  such  that  a  =  A  and  b  =  2A.  If  a 
and  b  are  increased,  the  computational  domain  in  an  FEM  simulation  increases.  This 
will  increase  the  number  of  unknowns  if  the  same  limits  on  MEL  are  used,  ultimately 
resulting  in  a  longer  solution  time.  The  MEL  can  be  increased  in  order  to  prevent  out- 
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Figure  5.6:  Larger  computational  domain 


of-memory  errors  during  solution  at  the  penalty  of  reduced  accuracy.  Increasing  the 
cloak  radii  in  the  Green’s  function  does  result  in  having  to  include  more  terms  in  the 
summation  due  to  the  constraints  of  Equation  5.2,  but  this  increase  in  computational 
budget  is  minimal  compared  to  the  increased  burden  in  an  FEM  simulation. 

As  an  example,  consider  the  simulation  geometry  shown  in  Figure  5.6.  More 
elements  are  going  to  be  needed  since  the  computational  domain  is  10A  x  10A.  A 
simulation  of  a  20-layer  cloak  of  homogeneous  material  approximating  the  material 
parameters  shown  Equation  5.16  was  performed  using  Comsol,  with  scattering  width 
results  obtained  and  compared  to  results  using  the  Green’s  function  formulation. 
These  results  are  shown  in  Figure  5.7.  The  A  between  the  two  simulations  was 
1.39  m2,  much  larger  than  the  0.1  m2  threshold.  As  before,  a  non-uniform  mesh  was 
applied  to  prevent  out-of-memory  errors.  The  MEL  for  the  layers  was  0.01A,  while  for 
the  remaining  areas,  MEL  =  0.3A.  Note  the  MEL  has  been  increased  compared  to  the 
same  20-layer  simulation  where  a  =  A  and  b  =  2A,  meaning  FEM  solution  accuracy 
will  decrease.  This  had  to  be  done  due  to  the  increase  in  the  size  of  the  computational 
domain.  The  resulting  mesh  consisted  of  972,698  elements  and  resulted  in  a  solution 
time  of  116  s.  In  addition  to  the  solution  time,  Comsol  took  168  s  to  simply  create 
the  mesh.  Further  increases  in  cloak  size  resulted  in  having  to  significantly  increase 
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FEM  and  GF  o2D  Comparison,  a  =  3A,,  b  =  4X 


Figure  5.7:  <T2d  for  larger  cylinder  and  cloak  size 

MEL  in  order  prevent  mesh  size  from  growing  beyond  the  computational  capabilities. 
The  increase  in  problem  geometry  had  little  impact  on  the  Green’s  function  solution. 
The  solution  time  did  increase  due  to  the  increase  in  the  number  of  terms  in  the 
summation,  but  the  solve  time  was  only  3.89  s.  The  Green’s  function  proved  to  be 
much  less  computationally  intensive  for  larger  problem  sizes 

Obviously  a  Green’s  function  approach  for  determining  scattering  widths  from 
a  cylindrical  cloak  results  in  a  significant  computational  savings.  The  computational 
domain  size  is  directly  related  to  the  cylindrical  cloaks’s  radius  in  that  a  larger  cloak 
results  in  a  larger  domain  size.  The  increase  in  computational  domain  requires  either  a 
longer  solution  time  due  to  the  increased  number  of  elements  or  a  reduction  in  mesh 
density  which  impacts  solution  accuracy.  The  Green’s  function  implementation  is 
much  faster  than  an  FEM  solution  and  is  more  adept  at  handling  problem  geometries 
which  require  denser  meshes  or  have  larger  computation  domains. 

5-4  Summary 

This  chapter  presented  a  Green’s  function  implementation  used  to  solve  for  the 
far  zone  scattered  held  from  a  two-dimensional  cylindrical  cloak.  The  Green’s  function 
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implementation  used  the  fact  thin,  concentric,  isotropic,  homogeneous  layers  could  be 
used  to  realize  a  cloak’s  anisotropic  properties.  The  Green’s  function  implementation 
was  shown  to  be  considerably  faster  when  solving  for  the  scattered  field  compared  to 
FEM  solutions  particulary  for  larger  computational  domains.  This  benefit  is  at  the 
expense  of  being  constrained  to  a  cylindrical  geometry.  A  different  Green’s  function 
would  have  to  be  derived  for  a  different  cloak  shape. 
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VI.  Metamaterial  Eigenfrequency  Decomposition 

Metamaterials  are  the  building  blocks  for  many  of  the  applications  defined  by  trans¬ 
formation  optics.  Transformation  optics  defines  a  material’s  constitutive  parameters 
necessary  to  achieve  a  desired  electromagnetic  field  behavior.  As  discussed  in  Chapter 
III,  the  unit  cells  which  make  up  a  metamaterial  are  designed  to  interact  with  the 
electromagnetic  fields,  resulting  in  a  desired  electric  or  magnetic  effect  (or  in  some 
cases  both).  The  unit  cell  metallizations  are  first  designed  with  knowledge  of  these 
field  interactions,  but  ultimately,  simulations  are  required  to  measure  the  exact  reso¬ 
nant  frequency  of  the  devices.  This  results  in  an  empirical  catalogue  of  measurements 
which  help  determine  how  changes  in  various  unit  cell  design  parameters  affect  the 
resonant  frequencies. 

Fischer  et  al.  used  an  eigendecomposition  to  design  substrates  for  patch  an¬ 
tennas  [31].  The  eigendecomposition  identifies  the  individual  eigenfrequencies  of  a 
structure.  These  eigenfrequencies  can  then  be  correlated  to  the  device  characteris¬ 
tics,  whether  it  be  geometry  or  material  property.  Fischer  et  al.  used  this  informa¬ 
tion  to  manipulate  eigenfrequency  location  by  changing  the  material  properties  of  the 
device.  This  resulted  in  an  increased  bandwidth  for  the  antennas.  A  similar  eigende¬ 
composition  design  method  could  be  applied  to  metamaterial  unit  cells.  This  initial 
investigation  determined  the  eigendecomposition  algorithm  is  applicable  to  metama¬ 
terial  unit  cell  designs.  What  follows  is  a  description  of  the  process  as  it  applies  to 
metamaterials  using  the  Cornsol  Multiphysics  software  package. 

6.1  Comsol  and  the  Finite  Element  Method 

The  Comsol  Multiphysics  software  package  is  used  to  obtain  all  FEM  solu¬ 
tions  in  this  chapter.  Therefore,  it  was  necessary  to  understand  exactly  how  Comsol 
implements  the  finite  element  method  to  solve  three  dimensional  electromagnetics 
problems.  Vector  basis  functions  (also  known  as  edge  elements)  are  used  within  each 
element  to  approximate  the  unknown  [27].  A  description  and  derivation  of  the  gen- 
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eral  form  for  the  vector  basis  functions  within  a  tetrahedral  subdomain  is  given  in 
Appendix  C. 

The  governing  partial  differential  equation  which  Cornsol  solves  is  the  vector 
wave  equation.  This  can  easily  be  derived  from  Maxwell’s  Equations.  Assuming  a 
source  free  domain,  Maxwell’s  equations  are 


V  x  E  =  —jojfiH, 


(6.1) 


V  x  H  =  jujeE. 


(6.2) 


The  electric  field  vector  wave  equation  can  be  developed  by  taking  the  curl  of  Equation 
6.1,  and  substituting  the  expression  for  V  x  H  from  Equation  6.2  into  this  expression. 
The  result  is 

-  )  —  k20erE  =  0.  (6.3) 

fir  J 

The  first  vector  Green’s  theorem  can  be  used  to  develop  the  weak  form  of  the  above 
equation  [95] 
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where  the  W i  terms  are  test  functions. 

The  FEM  approximates  E  within  each  geometric  element  using  vector  basis 
functions.  For  this  work,  tetrahedral  elements  are  used,  which  results  in  six  vector 
basis  functions  per  tetrahedral.  Within  each  element,  the  electric  field  is  approxi¬ 
mated  as 

6 

E  =  £  n]e‘.  (6.5) 

3= 1 

The  unknowns  are  the  Ee-  terms.  These  are  found  by  formulating  a  system  of  equations 
which  allows  the  unknowns  to  be  found.  The  system  of  equations  is  formed  by  first 
taking  the  dot-product  of  each  vector  basis  function  with  the  wave  equation,  and 
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integrating  the  result  over  the  volume  of  the  tetrahedral  element. 
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(6.6) 


Additionally,  the  test  functions,  W  i ,  are  defined  to  be  the  same  as  the  vector  basis 
functions  i.e  W i  =  iV*.  Finally,  replacing  E  in  Equation  6.6  with  the  approximation 
in  Equation  6.5  yields 


V  x  Ei  N,  ■  V  x  N, 
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(6.7) 


The  Ej’ s  are  constant  coefficients  and  can  be  removed  from  the  integration.  Addi¬ 
tionally,  the  dot-product  is  commutative;  thus,  the  order  can  be  reversed.  The  final 
equation  is 
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Note  that  each  integral  can  be  evaluated  analytically.  Additionally,  within  each  el¬ 
ement,  Hr  and  er  are  assumed  to  be  homogeneous  and  can  be  removed  from  the 
integration.  The  curl  terms  can  be  evaluated  based  on  the  following  [46]: 


Nx  =  WhiJtr 

V  x  Wili2  =  2VL?  x  X7Lei2. 


(6.9) 


The  If  represent  the  length  of  the  ith  edge  of  each  element  (Appendix  C  has  further 
details).  The  linear  interpolation  functions,  Lt] ,  Ll2 ,  i  =  1,2,  3, 4,  are  defined  as 


Ll  = 
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(6.10) 
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where  a®,  6®,  c®,  and  df  are  constants  defined  in  Appendix  C.  The  gradients  in  Equa¬ 
tion  6.9  are 

VL'  =  Tj  (6'x  +  c?$  +  dfz) ,  (6.11) 

e 

and  the  curl,  V  x  Ni,  can  be  written  as 
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Based  on  this  result,  it  is  possible  to  write  the  result  for  the  dot  product  of  the  two 
curls  as 
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The  above  are  all  constants.  Thus,  the  following  is  the  result  for  the  first  integral. 
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A  similar  procedure  can  be  done  to  determine  Ni  ■  Nj.  The  vector  basis  functions 
are  related  to  the  linear  interpolation  functions  by  [46] 


4  —  If  (LjjVI/jj  —  Li2VLtl) .  (6.15) 


^  e  e 

Therefore,  Ni  and  N  ■  can  be  expressed  as 

i  QVe  -^12  ^  d~  (-^21^22  -^22^1)  y  "f"  {Liidiz  -^22 ^21)  ? 

~^e  le 

N j  gye  [(-^ji^j2  —  x  +  {LjtCj2  —  Lj2Cj1 )  y  +  (KLj1dj2  —  Lj2 dj± )  z] . 

The  dot  product  can  be  expressed  by  the  following.  First  let 
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Then  the  dot  product  is 


/4e/e 

Ni  •  Nj  =  2  [LiiLi2fi2j2  ~  LhLj2fi2ji  —  Li2Lj1fi1j2  +  Li2Lj2fi1j1 \ 
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Note  the  L*  terms  are  functions  of  position.  Therefore,  the  evaluation  of  the  second 
integral  is  slightly  more  complicated.  The  following  formula  can  be  used  [46]. 
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As  an  example,  consider  Nl  ■  N1. 


e^2 


n[  ■  n[  =  -%-2  -  L!L2/21  -  L2L!/12  +  L2L2/n]  (6.20) 

(6Ve) 


Therefore,  the  second  integral  for  the  case  where  i  —  1  and  j  —  1  can  be  written  as 


y.  w  ■  Jvldv  =  [/22  /v,  Ui)  w  - 

/21  /ye  LiL2dV  —  f\2  JVe  L2LidV  +  /n  Jye  (L2)2dl/]  . 


(6.21) 


Using  the  integration  formula  and  noting  /12  =  /21,  the  result  is 


(6.22) 


It  is  possible  to  analytically  integrate  both  integrals  in  Equation  6.8.  Therefore, 
a  system  of  equations  can  be  formulated  such  that  the  unknown  terms,  UJ,  can  be 
determined.  Care  must  be  taken  when  constructing  the  global  matrices,  as  elements 
share  edges.  Knowing  the  method  used  to  construct  these  matrices  was  not  important 
to  this  effort.  Knowing  the  formulation  was  important  and  verified  as  such  in  the  next 
section. 
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6.2  Comsol  Formulation  Verification 

As  discussed  above,  the  FEM  converts  a  partial  differential  equation  into  a 
matrix  problem.  The  matrix  elements  consist  of  two  parts:  those  dependent  on  the 
permeability  and  those  dependent  on  permittivity.  Specifically,  examining  Equation 
6.8,  it  is  possible  to  define  the  following: 

A[f  =  —  f  V  x  JV-  •  V  x  N/V,  (6.23) 

Mr  Jve 

=  -k20er  [  Ni  ■  NjdV.  (6.24) 

Jve 

The  complete  impedance  matrix,  A,  is  defined  as 

A  =  A(1)  +  A(2).  (6.25) 

It  was  vital  to  be  able  to  access  the  matrices  defined  in  Equations  6.23  and  6.24.  How¬ 
ever,  in  its  default  settings,  Comsol  provides  access  to  only  the  complete  impedance 
matrix.  The  A ^  matrix  could  be  obtained  by  setting  all  domain  permittivities  to 
zero,  effectively  eliminating  the  A ^  contribution.  The  A ^  matrix  could  then  be 
found  by  first  solving  for  A  and  then  subtracting  the  A ^  matrix. 

The  kb2)  matrix  could  also  be  directly  evaluated  by  manually  modifying  the 
partial  differential  equation  which  Comsol  uses  to  generate  the  system  of  equations 
[34],  Effectively,  the  governing  weak  form  of  the  partial  differential  equation  was 
changed  to 

—A;2  /  erE  ■  N fiV  =  0.  (6.26) 

Jve 

Note  Equation  6.26  is  scaled  by  k2.  Therefore,  checking  the  validity  of  the  A ^  matrix 
was  possible.  Doubling  kQ  resulted  in  each  term  in  the  A )2)  matrix  being  scaled  by  a 
factor  of  four.  Additionally,  the  kb2)  matrix  obtained  by  altering  the  equation  system 
was  compared  to  that  obtained  by  simply  solving  for  the  difference  of  A  and  kb1). 
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The  results  were  the  same  to  within  10  10  for  various  test  cases.  Therefore,  it  was 
concluded  the  correct  A ^  and  A ^  matrices  could  be  obtained  from  Comsol. 

6. 3  Eigendecompostion 

A  standard  FEM  formulation  will  solve  a  system  of  equations  of  the  form 

+  A<2)(fc)][£7(fc)]  =  [/(fc)],  (6.27) 

where  [E(k)\  are  the  unknowns,  [f(k)\  is  some  forcing  function,  and  k  is  the  wave 
number  of  the  excitation.  Note  wave  number  is  directly  related  to  frequency. 

The  impedance  matrix,  A,  has  two  components,  A^  and  A^2\  Note  from 
Equations  6.23  and  6.24  that  A(i>  is  independent  of  frequency  ( k ),  while  A ^  is 
dependent  on  it.  Hence,  a  normalized  wave  number,  k ,  can  be  defined  as  [32] 

k  =  (6.28) 

kQ 

Equation  6.28  can  be  used  to  write  an  equivalent  expression  for  Equation  6.27. 

+  kAW(k„)}{E(k)]  =  [/(*)]  (6.29) 

The  unknowns  in  Equation  6.29  can  be  found  by 

[E{k)]  =  [A^(k0)  +  M2\k0)}^[f(k)}.  (6.30) 

The  normalized  eigenfrequencies  are  those  values  of  k  which  make  the  matrix  inversion 
on  the  right  side  of  Equation  6.30  impossible.  The  eigenfrequencies  can  be  analytically 
found  using  the  method  put  forth  by  Fischer  et  al.  First  an  eigendecomposition  on 
the  A^  and  A matrices  is  performed. 

XAX(~V  =  [A(2)(/c0)](-1)A(1)(/c0)  (6.31) 
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X  is  a  matrix  whose  columns  are  the  eigenvectors  while  A  is  a  diagonal  matrix  whose 
elements  are  the  eigenvalues.  An  equivalent  expression  for  Equation  6.30  can  be 
written  using  the  eigendecomposition. 


[E{k)} 


k2I  +  X  AX(”1} 


(-1) 


[A(2,(fc0)]<-1)  [/(*)] 


(6.32) 


The  I  term  is  the  identity  matrix.  Equation  6.32  can  be  rewritten  by  defining 

A~k  =  A  +  D(k2),  (6.33) 

where  D{k2)  is  a  diagonal  matrix  with  elements  k2.  Therefore,  an  equivalent  form  for 
Equation  6.30  is  [32] 


[*W]  =  XA<-1)X<-1>p<2)(*o)](-1)[/(*:)]-  (6.34) 

Note  that  AX)  1W^-1*[A*2^(fc0)]'“6  is  resonant  at  kj  —  —Re{ A,).  This  is  because 
the  elements  of  the  diagonal  matrix  A~k  are  A*  +  k2.  Therefore,  by  defining  ki  as  such, 
there  becomes  a  zero  on  the  diagonal  of  the  Ak  matrix,  making  it  noninvertible  or 
resonant  [31].  Thus,  the  eigenfrequencies  for  Equation  6.29  can  be  determined. 


6.4  Eigendecomposition  Verification 

In  order  to  verify  the  Comsol-generated  A ^  and  A ^  matrices  could  be  used  to 
identify  the  eigenfrequencies  for  structures  by  the  method  described  in  the  previous 
section,  a  simple  test  case  was  performed.  A  rectangular  cavity  as  shown  in  Figure 
6.1  has  easily  determined  eigenvalues.  When  PEC  boundary  conditions  are  applied 
to  all  six  faces,  the  TE~  eigenvalues  are  given  by  [13] 


(fr) 


mnp 


(6.35) 
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Figure  6.1:  PEC  Rectangular  Resonator  Geometry 


where  m  —  0,1,  2...,  n  =  0, 1,  2...,  p  —  1,  2,  3...  and  m  =  n  ^  0.  The  TM"  modes  are 
given  by  the  same  equation  with  the  exception  that  m  =  1,2,3...,  n  =  1,2,3...,  and 
=  0,1,2.... 

A  PEC  rectangular  resonator  as  shown  in  Figure  6.1  was  created  in  Cornsol 
where  a  =  0.02286  m,  b  =  0.01016  m,  and  c  =  0.1143  m.  The  eigenfrequencies 
were  extracted  using  the  method  described  in  the  previous  section.  These  are  plotted 
in  Figure  6.2.  The  top  graph  shows  all  extracted  eigenfrequencies  while  the  bot¬ 
tom  graph  zooms  in  to  the  region  where  the  eigenfrequencies  from  6-20  GHz  are 
shown.  For  the  rectangular  resonator  shown  in  Figure  6.1  with  values  previously  de¬ 
fined  for  a,  b,  and  c,  the  theoretically  calculated  eigenfrequencies  were  compared  to 
the  extracted  eigenfrequencies.  These  results  are  shown  in  Table  6.1.  The  extracted 
eigenfrequencies  matched  the  theoretical  values  to  within  0.25%  for  all  21  eigenfre¬ 
quencies  less  than  16  GHz.  Above  16  GHz,  the  extracted  frequencies  did  match  the 
theoretical  ones.  However,  the  extraction  also  produced  non-theoretical  values.  For 
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Extracted  Eigenfrequencies 
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Figure  6.2:  PEC  Resonator  Eigenfrequencies 

example,  there  are  50  theoretical  eigenfrequencies  below  20  GHz.  The  extraction  al¬ 
gorithm  produced  61.  There  are  145  eigenfrequencies  below  30  GHz.  The  extraction 
algorithm  produced  205. 

Equation  6.35  shows  there  to  be  an  infinite  number  of  eigenfrequencies.  The 
number  of  extracted  eigenfrequencies  is  finite  because  the  number  of  eigenfrequencies 
is  related  to  the  size  of  the  impedance  matrices.  It  was  found  the  frequency  at 
which  the  algorithm  begins  to  extract  non-theoretical  values  is  also  related  to  the 
size  of  the  impedance  matrix.  The  size  of  the  impedance  matrices  can  be  controlled 
by  increasing  the  mesh  fidelity  in  Cornsol.  However,  increases  in  mesh  fidelity  will 
result  in  significant  increases  in  computation  time  because  the  eigendecompostion 
shown  in  Equation  6.31  is  an  0(N3)  operation  [32],  As  an  example,  the  results 
shown  in  Figure  6.2  were  obtained  with  a  9,516  x  9,516  element  impedance  matrix. 
Eigenfrequency  extraction  execution  time  was  75  minutes.  A  less  dense  mesh  resulting 
in  an  impedance  matrix  with  1,190  x  1,190  elements  resulted  in  a  solution  time  of  7.57 
seconds.  However,  the  extraction  algorithm  found  nontheoretical  eigenfrequencies  at 
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Table  6.1:  Eigenfrequencies  in  GHz  for  PEC  Rectangular  Resonator 


m 

n 

P 

f Theory 

f Extracted 

%  Error 

1 

0 

1 

6.687 

6.686 

0.009 

1 

0 

2 

7.062 

7.061 

0.018 

1 

0 

3 

7.647 

7.646 

0.009 

1 

0 

4 

8.397 

8.394 

0.030 

1 

0 

5 

9.273 

9.268 

0.050 

1 

0 

6 

10.243 

10.237 

0.049 

1 

0 

7 

11.281 

11.273 

0.076 

1 

0 

8 

12.372 

12.362 

0.085 

2 

0 

1 

13.180 

13.169 

0.082 

2 

0 

2 

13.374 

13.363 

0.084 

1 

0 

9 

13.502 

13.481 

0.152 

2 

0 

3 

13.692 

13.678 

0.099 

2 

0 

4 

14.125 

14.105 

0.136 

1 

0 

10 

14.662 

14.638 

0.161 

2 

0 

5 

14.662 

14.646 

0.110 

0 

1 

1 

14.812 

14.774 

0.253 

0 

1 

2 

14.985 

14.949 

0.235 

0 

1 

3 

15.269 

15.240 

0.189 

2 

0 

6 

15.294 

15.275 

0.120 

0 

1 

4 

15.659 

15.627 

0.202 

1 

0 

11 

15.846 

15.823 

0.143 

2 

0 

7 

16.008 

15.985 

0.145 

12.79  and  13.00  GHz.  Additionally,  the  error  between  the  extracted  and  theoretical 
eigenfrequencies  increases  as  mesh  fidelity  decreases.  For  example,  the  error  for  the 
101  mode  for  the  reduced  mesh  is  0.33%,  still  small,  but  much  larger  than  the  result  in 
Table  6.1.  Also,  the  extraction  for  the  less  dense  mesh  produced  72  eigenfrequencies 
less  than  20  GHz.  Obviously  better  results  are  obtained  with  a  denser  mesh  at  the 
penalty  of  increased  solve  times. 


6.5  S-Parameter  Measurements 

As  discussed  in  Chapter  III,  the  published  literature  shows  no  S-parameter 
results  had  been  obtained  using  the  Cornsol  Multiphysics  software.  To  validate  the 
software’s  capability,  a  unit  cell  as  described  in  [90]  was  created  in  Cornsol.  This  unit 
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cell  is  shown  in  Figure  6.3.  The  dimensions  of  the  unit  cell  were  identical  to  those 


Figure  6.3:  Unit  Cell  for  S-Parameter  Measurements 


used  by  Smith  et  al.  The  cell  size  was  cubic  with  d  =  2.5  mm.  The  substrate  was 
0.25  mm  thick  with  er  =  4.4  —  j0.088.  The  outer  ring  was  2.2  mm  with  a  linewidth 
of  0.2  mm.  The  gap  between  the  inner  and  outer  ring  was  0.15  mm,  and  the  gap 
between  the  ring  end’s  were  each  0.3  mm.  The  rod  on  the  opposite  side  had  a  width 
of  0.14  mm  and  was  2.5  mm  long.  All  metal  thicknesses  were  17  pm  and  were  given 
the  properties  of  a  PEC.  S-parameter  measurements  were  simulated  in  Comsol.  The 
constitutive  parameters  were  extracted  using  the  method  described  in  Section  3.3.2. 
All  results  were  compared  to  the  published  ones  and  are  shown  in  Figure  6.4.  The 
Comsol ’s  results  were  nearly  identical  to  those  in  [90].  Therefore,  it  was  concluded 
the  software  could  accurately  provide  S-parameter  measurements. 

The  goal  of  this  chapter  is  to  perform  an  eigenfrequency  extraction  on  the  unit 
cell  shown  in  Figure  6.3.  However,  the  size  of  the  impedance  matrix  would  be  an 
issue.  The  impedance  matrix  resulting  from  a  Comsol  S-parameter  simulation  of 
the  unit  cell  shown  in  Figure  6.3  was  352,692  x  352,692.  Using  even  the  coarsest 
settings  to  create  the  mesh  resulted  in  a  39,142  x  39,142  impedance  matrix.  Various 
tests  with  the  eigendecomposition  algorithm  showed  any  matrices  larger  than  20,000 
x  20,000  elements  would  overwhelm  the  machine.  Therefore,  the  thickness  of  the 
metal  was  eliminated  by  making  all  metal  structures  infinitely  thin  PEC  boundaries. 
There  was  a  concern  this  would  drastically  change  the  resonant  characteristics  of 
the  structure.  However,  S-parameter  measurements  were  taken  and  compared  to  the 
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Figure  6.4:  S-Parameter  Measurements  and  Extracted  Index  of  Refraction 


Table  6.2:  Mesh  Density  and  Impedance  Matrix  Size 


Mesh  Setting 

Mesh  Number 

Impedance  Matrix  Size 

Normal 

One 

116,932  x  116,932 

Coarse 

Two 

34,578  x  34,578 

Coarser 

Three 

20,866  x  20,866 

Extra  Coarse 

Four 

7,248  x  7,248 

Extremely  Coarse 

Five 

3,260  x  3,260 

original  results.  Infinitely  thin  metallic  structures  slightly  shifted  the  S-parameter 
measurements,  but  by  approximately  0.1  GHz.  There  was  still  the  same  resonant 
type  behavior  when  using  an  embedded  PEC  structure  compared  to  PEC  structures 
with  small  thicknesses.  The  advantage  of  removing  the  metal  thickness  is  significantly 
fewer  mesh  elements  are  required. 

Using  embedded  PEC  structures  resulted  in  a  smaller  impedance  matrix.  How¬ 
ever,  the  matrices  were  still  too  large.  The  density  of  the  mesh  was  varied  to  determine 
mesh  density’s  impact  on  the  S-parameter  measurements.  Five  different  mesh  set¬ 
tings  were  used,  each  resulting  in  different  impedance  matrix  sizes.  These  are  shown 
in  Table  6.2  with  the  resulting  S-parameter  measurements  shown  in  Figure  6.5. 
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S- Parameter  Comparison 


Figure  6.5:  S-Parameter  Magnitude  Comparison  for  Mesh  Densities 

Based  on  these  results,  it  is  obvious  the  impact  of  using  a  smaller  impedance 
matrix  (less  dense  mesh)  does  not  significantly  impact  the  resonant  characteristics 
of  the  unit  cell.  In  the  previous  section  it  was  shown  the  size  of  the  impedance 
matrix  directly  relates  to  the  accuracy  of  the  extracted  eigenfrequencies.  However, 
due  to  computer  limitations,  the  densest  impedance  matrix  used  for  the  unit  cell 
eigendecomposition  in  the  following  section  was  mesh  number  four,  which  had  an 
impedance  matrix  of  7,248  x  7,248. 

6.6  Unit  Cell  Eigendecomposition 

An  eigendecomposition  was  performed  on  unit  cells  having  the  ring  characteris¬ 
tics  shown  in  Figure  6.3.  To  help  further  reduce  the  impedance  matrix  size,  the  metal 
rod  was  removed  from  the  structure.  This  rod  impacts  the  effective  permittivity  of 
the  unit  cell  structure.  The  permeability  is  a  function  of  the  metal  rings. 

The  ring  characteristics  were  slightly  modified.  The  gap  between  the  ring  ends 
was  changed  to  0.1  mm  (Mod  1)  and  0.6  mm  (Mod  2).  This  significantly  changes  the 
capacitance  of  the  unit  cell  structure,  which  will  alter  the  resonant  behavior.  This  can 
be  seen  in  the  S-parameter  measurements  shown  in  Figure  6.6.  Note  the  pronounced 
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Gap  Modification  S-Parameters 


Figure  6.6:  S-Parameter  Magnitudes  for  Gap  Modifications 

shift  in  the  S-parameters  by  more  than  1  GHz.  An  eigendecompostion  was  done  on  the 
Mod  1  and  Mod  2  configurations.  The  eigenvalues  for  each  geometry  were  extracted 
using  the  eigendecompostion  method  described  in  this  chapter.  Plots  of  the  extracted 
eigenvalues  from  each  structure  and  the  magnitude  of  their  differences  are  shown  in 
Figure  6.7.  Note  each  configuration  resulted  in  128  nonzero  eigenvalues  less  than  100 
GHz.  Additionally,  some  values  were  significantly  different,  but  these  large  differences 
manifested  in  the  higher  eigenfrequencies.  Those  eigenvalues  in  the  6-18  GHz  range 
were  very  close.  This  is  better  seen  in  Figure  6.8.  Based  on  the  S-parameter  measure¬ 
ments,  the  structures  have  significantly  different  resonant  behaviors.  However,  this  is 
not  as  obvious  when  examining  the  individual  eigenfrequencies.  It  does  appear  there 
are  specific  resonant  regions.  Note  the  lack  of  eigenvalues  between  5.9  and  6.3  GHz, 
8.5  and  9.4  GHz,  and  13.4  and  14.1  GHz.  However,  no  definitive  conclusions  could  be 
drawn  from  this  data.  Further  eigenfrequency  extractions  were  performed  with  the 
gap  in  the  rings  modified  to  different  widths  (0.2,  0.3,  0.4,  and  0.5  mm).  Extracted 
eigenfrequencies  were  compared,  but  no  definitive  changes  in  specific  eigenfrequencies 
were  noted. 
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Extracted  Eigenvalues 


Eigennumber 


Figure  6.7:  Eigenvalues  for  Gap  Modifications 


It  was  determined  the  reduced  impedance  matrix  sizes  were  not  allowing  enough 
fidelity  to  see  the  eigenvalue  changes.  To  test  this,  two  different  meshes  were  used  to 
create  different,  although  similarly  sized,  impedance  matrices  for  the  Mod  2  structure. 
The  different  meshes  were  created  by  toggling  the  jiggle  function  within  Comsol.  The 
extracted  eigenfrequencies  obtained  for  the  same  structure  but  using  different  meshes 
were  compared.  Results  are  shown  in  Figure  6.9.  The  extracted  eigenfrequencies 
are  not  identical.  In  fact,  note  the  similarities  between  Figures  6.7  and  6.9.  The 
graphs  of  the  differences  between  the  extracted  eigenfrequencies  in  both  figures  seem 
to  indicate  the  variations  in  the  eigenfrequencies  are  a  function  of  the  mesh  rather 
than  a  function  of  changes  in  the  structure  of  the  unit  cell.  Unfortunately,  increasing 
the  mesh  size  to  the  maximum  capability  did  not  alleviate  this  problem.  It  is  believed, 
however,  that  further  increases  in  mesh  fidelity  will  lead  to  specific  eigenfrequency 
identification  of  different  unit  cell  structures.  However,  this  could  not  be  tested  due 
to  memory  limitations  on  the  computational  hardware. 
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Extracted  Eigenvalues 
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Figure  6.8:  Eigenvalues  from  6-18  GHz 

6. 7  Summary 

This  chapter  presented  an  eigenfrequency  decomposition  method  which  has  been 
shown  to  be  able  to  extract  a  structure’s  individual  eigenfrequencies.  The  method 
was  implemented  using  the  Comsol  Multiphysics  software  and  tested  using  a  structure 
with  theoretically  known  eigenfrequencies.  Extracted  eigenfrequencies  of  a  rectangu¬ 
lar  PEC  resonator  matched  the  theoretical  values.  However,  it  was  shown  mesh 
density  plays  an  important  roll  in  the  fidelity  of  the  solution. 

The  eigenfrequency  decomposition  method  was  applied  to  a  metamaterial  unit 
cell.  In  order  to  reduce  impedance  matrix  size,  infinitely  thin  metallic  boundaries 
were  used  in  place  of  actual  metal  structures.  S-parameter  measurements  showed 
this  change  had  little  impact  on  the  device’s  resonant  behavior.  Eigenfrequencies 
were  extracted  from  the  unit  cells.  However,  the  mesh  was  not  dense  enough  to 
allow  identification  of  shifts  in  eigenfrequencies  as  a  result  of  changes  to  the  device 
structures.  The  extracted  eigenfrequencies  were  impacted  by  changes  in  the  mesh. 
A  finer,  denser  mesh  is  needed  to  adequately  simulate  these  unit  cells.  However, 
computational  limits  were  reached  due  to  memory  limitations. 


Gap  =  0.1  mm 
Gap  =  0.6  mm 


Eigennumber 


111 


100 


|  50 
o 
n 
cr 
2 
LL 

0 

0  20  40  60  80  100  120 

Eigennumber 


Extracted  Eigenvalues 


Mesh  1 


100 


50 


0 


Mesh  2 


X 


0  20  40  60  80  100  120 


Eigennumber 


Figure  6.9: 


Extracted  Eigenfrequencies  Using  Different  Meshes 
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VII.  Conclusions 


7.1  Research  Summary 

The  performance  of  simplified  cylindrical  cloaks  with  various  constitutive  pa¬ 
rameters  was  analyzed  in  order  to  understand  the  impact  constitutive  values  have  on 
field  behavior.  Prior  to  this  research,  the  material  parameters  of  simplified  cloaks  have 
focused  on  satisfying  specific  values  of  £zyg  and  ezyr  while  matching  the  impedance 
at  the  cloak’s  outer  boundary.  A  third  constraint  equation  was  introduced  which 
helps  control  the  overall  effectiveness  of  the  cylindrical  cloak. 

Cylindrical  cloaks  were  analyzed  with  constitutive  parameters  that  satisfied 
the  specific  values  for  £zyg  and  ezyr.  It  was  shown  deviations  from  this  derived 
third  constraint  equation  resulted  in  larger  fields  being  transmitted  into  a  cylindrical 
cloak’s  hidden  region.  As  the  cloak’s  constitutive  parameters  were  changed  such  that 
this  new  constraint  was  better  satisfied,  the  amount  of  energy  transmitted  into  the 
hidden  region  was  shown  to  be  reduced.  The  resulting  impedance  mismatch  at  r  —  b 
due  to  changing  the  constitutive  parameters  resulted  in  a  significant  scattered  field. 
However,  despite  reducing  energy  transmitted  into  the  hidden  region,  which  resulted 
in  a  reduction  in  the  scattered  field  by  the  cloaked  object,  the  cloak  itself  was  creating 
a  large  scattered  field.  Hence,  in  terms  of  overall  scattering  width,  having  a  matched 
impedance  at  r  —  b  was  shown  to  be  more  important  than  reducing  the  transmitted 
energy  into  the  hidden  region. 

A  new  way  to  develop  simplified  material  parameter  sets  for  cylindrical  cloaks 
was  developed.  Specifically,  for  TMZ  incident  waves,  the  approximation  of  yg  should 
first  be  defined  using  a  Taylor  series  expansion  of  the  ideal  parameter  as  defined  by  the 
derived  third  constraint  equation.  The  constitutive  parameters  yr  and  can  then  be 
determined  by  making  the  products  yg£z  and  yr£z  equal  to  the  same  products  using 
the  ideal  material  parameter  set.  The  performance  of  cloaks  developed  in  this  manner 
is  limited  only  by  the  number  of  terms  used  in  the  Taylor  series  expansion,  which  is 
dictated  by  existing  manufacturing  capabilities.  Additionally,  the  applicability  of  this 
method  extends  to  TES  fields  by  duality. 
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Scattering  width  improvement  was  observed  for  all  angles  when  compared  to 
previous  published  material  parameter  sets.  Significant  improvement  was  noted  in 
the  forward  scattering  region.  It  was  also  shown  the  simplified  parameter  set  put 
forth  in  [102]  is  a  simplification  of  this  method  in  which  the  Taylor  series  expansion 
of  no  is  limited  to  the  first  term.  These  parameter  sets  were  found  to  have  relatively 
consistent  performance  for  all  values  of  the  cloak’s  radial  boundary,  b.  Performance 
for  a  constant  number  of  terms  in  the  Taylor  series  does  slightly  degrade  as  b  increases, 
but  for  all  6,  ideal  cloaking  performance  is  approached  as  IV  — >  oo. 

A  Green’s  function  approach  for  determining  scattering  widths  from  a  cylindri¬ 
cal  cloak  was  shown  to  have  significant  computational  savings  compared  to  standard 
FEM  methods.  This  savings  can  be  useful  for  error  analysis  or  optimization  studies  on 
a  particular  cloak  geometry.  Also,  the  computational  domain  size  is  directly  related 
to  the  cylindrical  cloak’s  radius.  A  larger  cloak  results  in  a  larger  domain  size.  The 
increase  in  computational  domain  requires  either  a  longer  solution  time  due  to  the 
increased  number  of  elements  or  a  reduction  in  mesh  density  which  impacts  solution 
accuracy.  The  Green’s  function  implementation  is  much  faster  than  an  FEM  solution 
and  is  more  adept  at  handling  larger  problem  geometries. 

Metamaterial  unit  cells  were  analyzed  using  an  eigendecomposition  technique. 
S-parameter  measurements  showed  definite  shifts  in  unit  cell  resonant  frequencies  due 
to  structural  changes.  Eigenfrequencies  were  extracted  from  the  unit  cells.  Shifts  in 
resonant  frequency  locations  were  noted  for  different  cell  geometries,  but  no  defini¬ 
tive  relationships  could  be  drawn.  Mesh  densities  were  limited  to  very  coarse  set¬ 
tings  due  to  computational  limitations.  Extracted  eigenfrequencies  were  shown  to  be 
mesh-dependent.  The  problem  could  be  ameliorated  by  increasing  mesh  fidelity,  but 
memory  limitations  were  reached  preventing  this  was  being  further  explored. 

7.2  Recommendations  for  Future  Research 

The  work  done  in  Chapter  IV  considered  infinite  cylindrical  cloaks.  Such  anal¬ 
ysis  was  done  due  to  computational  efficiencies  gained  when  solving  two-dimensional 
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problems.  To  physically  realize  such  a  structure,  a  three-dimensional  analysis  will  be 
required  i.e.  the  cloak  must  be  terminated  in  the  5  direction  at  z  =  Z\  and  z  =  z2. 
This  is  shown  in  Figure  7.1.  The  termination  of  the  infinite  cylinder  will  result  in 

*  z 


^ Z  =  z: 

Figure  7.1:  Three-Dimensional  Cylindrical  Cloak 


scattering  from  the  cloak  even  if  the  ideal  parameters  are  assumed.  This  is  due  to 
trailing  edge  diffraction  resulting  from  the  edges  at  Z\  and  z 2  in  Figure  7.1.  ft  is 
not  clear  how  large  the  diffracted  field  will  be.  The  typical  2-D-to-3-D  conversion 
formula  [49], 

2 12 

&3D  =  a2  D~^i  (7-1) 

cannot  be  used  since  a2D  for  ideal  cloaks  is  zero.  However,  for  reduced  parameter 
cloaks,  it  will  be  interesting  to  analyze  the  accuracy  of  Equation  7.1.  The  z  directed 
incident  held  will  excite  surface  currents  in  the  same  direction.  There  are  sharp 
discontinuities  at  z  —  Z\  and  z  =  z2,  which  will  result  in  a  significant  scattered  held. 
The  size  of  the  scattered  held  will  depend  on  the  properties  of  the  terminating  ends 
of  the  cloak.  Ways  to  reduce  to  the  scattering  could  include  tapering  using  a  resistive 
material  to  better  match  the  termination  to  free  space. 
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Additionally,  a  three-dimensional  cloak’s  functionality  with  the  ideal  cloak  pa¬ 
rameters  should  be  independent  of  9.  As  the  incident  angle  is  swept  around  the  body 
of  the  cloak.  The  diffraction  effects  should  be  identical  for  all  incident  6  angles.  How¬ 
ever,  when  9  is  changed  i.e.  the  incident  wave  approaches  from  a  direction  not  normal 
to  the  cylinder’s  broad-side,  there  will  be  a  significant  change  in  RCS  because  the 
terminated  ends’  scattering  mechanism  is  changing  from  that  of  a  grazing  incidence 
to  a  specular  return. 

This  work  focused  specifically  on  two-dimensional  cylindrical  cloaks.  A  compu¬ 
tational  improvement  was  noted  by  using  a  Green’s  function  to  compute  the  scattering 
width  of  a  cloaked  PEC  cylinder.  A  PEC  has  a  complex  permittivity  such  that 

£c(r,(/),u )  =  lim  £  +  j—„  (7.2) 

(7  — »00  (jJ 

where  o  is  the  conductivity  of  the  material.  Alternate  Green’s  functions  can  be 
derived  such  that  o  is  finite  at  r  =  a,  which  results  in  field  penetration  into  the 
hidden  region.  The  material  inside  the  hidden  region  may  be  inhomogeneous  and 
not  symmetric  with  respect  to  9.  This  would  result  in  a  significantly  more  complex 
Green’s  function,  but  allows  for  any  geometry  to  be  placed  inside  the  inner  boundary. 

As  discussed  in  Chapter  II,  there  have  been  a  number  of  different  cloak  geome¬ 
tries  discussed  in  the  published  literature.  It  would  be  an  interesting  academic  exercise 
to  derive  the  Green’s  function  for  these  geometries  and  implement  a  computational 
solution  as  was  done  in  Chapter  V.  Green’s  functions  can  help  gain  physical  insight 
into  the  cloaking  function,  which  may  prove  useful  for  these  alternate  geometries.  Ad¬ 
ditionally,  it  is  expected  the  Green’s  function  will  provide  a  significant  computational 
improvement  which  may  be  beneficial  if  optimizations  are  being  performed. 

Finally,  further  work  can  be  done  extracting  the  eigenfrequencies  from  unit  cell 
designs.  Specifically,  increasing  the  memory  should  allow  for  denser  meshes.  This 
increase  in  mesh  density  should  allow  for  structural  differences  to  manifest  in  eigen- 
frequency  shifts.  This  could  create  a  new  unit  cell  design  paradigm  where  eigenfre- 
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quency  location  is  correlated  to  a  structural  feature.  Designs  could  be  optimized  for 
larger  bandwidths  by  using  Fischer’s  optimization  technique  described  in  [31].  Addi¬ 
tionally,  alternate  techniques  to  address  the  bandwidth  problem  of  unit  cells  should 
be  investigated.  Specifically,  multiresonant  structures  within  the  same  unit  cell  might 
provide  an  increase  in  bandwidth.  Additionally,  active  materials  could  be  used  in  unit 
cell  designs  to  increase  or  actively  change  the  bandwidth  characteristics  of  the  unit 
cells. 
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Appendix  A.  Material  Parameter  Derivation 

The  development  below  mirrors  what  is  shown  in  [97]  with  clarifying  text  added  where 
the  original  text  was  deemed  ambiguous  or  unclear.  Maxwell’s  equations  govern  the 
behavior  of  electromagnetic  waves;  they  are  valid  in  any  coordinate  system.  Ward 
and  Pendry  have  shown  the  behavior  of  electromagnetic  waves  in  a  general  coordinate 
system  can  be  modeled  in  Cartesian  coordinates  using  a  designed  material  with  spe¬ 
cific  permittivity  and  permeability  tensors  [72,97].  It  is  this  fact  upon  which  cloaking 
is  based.  To  prove  this,  take  the  general  form  of  Faraday’s  Law  in  free  space. 

VxE  =  - ii0 —  (A.l) 

The  desire  is  to  find  the  form  of  Faraday’s  law  in  a  general  coordinate  system  given  by 
the  variables  (qi,  q2,  q3)  with  unit  vectors,  U\,  u2,  and  u3  in  the  direction  of  the  g1;  q2, 
and  q3  axes.  Additionally,  it  is  assumed  there  exists  a  transformation  from  Cartesian 
coordinates  to  this  general  coordinate  system  where  the  point  q%  is  a  function  of 
(x,  y,  z )  expressed  as 

Qi  =  F1  (x,y,z), 

q2  =  F2(x,y,z),  (A.2) 

<?3  =  F3(x,y,z). 

It  is  assumed  the  transformation  is  invertible  such  that  the  point  (<q ,  q-.i)  can  be 

transformed  back  to  Cartesian  coordinates  by 

x  =  /l(gi,<?2,<?3), 

y  =  /2  (51,92,93),  (A-3) 

Z  =  Mqi,Q2,q3)- 

Next,  it  is  important  to  understand  how  to  calculate  the  differential  length,  ds,  of  a 
line  segment  in  the  general  coordinate  system.  This  requires  the  use  of  the  Euclidean 
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metric,  which  is  defined  in  terms  of  the  Jacobian,  J,  and  is  equal  to  JTJ  where 
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(A.4) 


The  squared  length  of  a  differential  line  segment  is  expressed  as 


ds 2 


3  3 


(A.5) 


Obviously,  in  Cartesian  coordinates,  q\  =  x,  q2  =  y,  and  q3  =  z  and  the  resulting 
Jacobian,  transpose  Jacobian,  and  product  of  the  two  matrices  are 


/  i  n  n  \ 


J  = 
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1  0  0 
0  1  0 
0  0  1 
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Hence,  the  squared  differential  length  is  the  familiar  form 


V 


/ 


(A. 6) 


(A.7) 


ds2  =  dx2  +  dy 2  +  dz2. 


(A.8) 


In  their  paper,  Ward  and  Pendry  use  the  following  notation 


ds 2  —  Qndq2  +  Q22d^2  T  Qz^dq^  +  2Q \2dqidq-2  +  2Q i^dqidq-^  +  2Q2^dq2dq 3,  (A. 9) 


where 


Qij 


dx  dx  dy  dy  dz  dz 

dqt  dqj  dq{  dqj  dqt  dqj  ' 


(A.10) 


119 


Using  the  definition  of  Qt]1  it  is  obvious  Equation  A. 9  matches  that  of  Equation  A. 5. 

It  is  necessary  to  represent  the  differential  length  of  a  line  segment  in  the  di¬ 
rection  of  one  of  the  unit  vectors.  As  an  example,  suppose  there  exists  a  differential 
length  only  in  the  U\  direction.  The  result  is  dq-2  =  dq-z  =  0.  This  results  in  the 
length,  ds  to  be 

ds  =  \J  (JT  J)11dqi.  (A.  11) 

Using  the  notation  of  Ward  and  Pendry,  this  is  identical  to 

ds  =  y/Qndqi.  (A.  12) 

As  a  simplification,  Ward  and  Pendry  let  Q\  =  Q n  and  Q i  =  y/Qn.  Therefore 

ds  =  Q\dq\,  (A. 13) 

or,  in  more  general  terms,  the  differential  length  along  the  direction  of  the  ith  unit 
vector  is 

dsi  =  Qidqt .  (A.  14) 

This  notation  is  necessary  because  the  desire  is  to  determine  the  form  of  Faraday’s 
Law  in  the  general  coordinate  system.  To  do  this,  consider  a  small  differential  element 
of  the  shape  of  a  parallelepiped  as  shown  in  Figure  A.l.  It  is  possible  to  calculate 


Figure  A.l:  Differential  parallelepiped  element 
the  projection  of  V  x  E  onto  the  normal  to  the  ■u1-u2  plane.  To  do  this,  first  take  the 
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line  integral. 


E-dl 


(A.  15) 


c 


The  integral  can  be  evaluated  by  letting  Et  =  E  ■  ut .  The  components  along  the  closed 
path  are  shown  in  Figure  A. 2.  By  taking  into  account  the  direction  of  the  contour, 


A  If/A; 


EiQidch + —[EiQidlh\  dch 
dq2 

— - 


E2i 22^2  +  ^  [E2Q2dq2\dq l 


EiQA i 

Figure  A. 2:  Line  integral  differential  components 


the  line  integral  is  evaluated  to  be 


-*  -a-  Q  Q 

E-dl  =dqi  —  [ E2Q2dq2 }  ~  [E^Q^qi]  . 

C  vqi  c>q2 


(A. 16) 


Stokes  Theorem  can  now  be  applied.  Recall,  Stokes  Theorem  states 


E  ■  dl  =  /  /  V  x  E  ■  hdS  =  /  /  V  x  E  ■  (tR  x  u2 )  QidqiQ2dq2-  (A. 17) 

C  J  J  S  J  J  s 


For  this  geometry  note  that 


n  —  u i  x  u2, 
dS  =  QidqiQ2dq2. 


(A.  18) 


Thus,  applying  Stokes  Theorem  to  the  line  integral  along  this  differential  contour 
yields 


-t  d  d 

V  x  E  •  (hi  x  u2)  Q!dqiQ2dq2  =  dq1—~  [ E2Q2dq2 \  -  dq2——  [FdQidgi]  •  (A.  19) 

oqi  oq2 
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This  can  be  simplified  to 


o  o 

v  x  e  •  («!  x  fi2)  g^  =  —  [E2g2]  -  —  [e1q1]  .  (A.20) 

dq2 

Letting  E{  =  QiEi  results  in  the  following. 

LL1  r)  fi1 

V  x  £  •  (ux  x  u2)  QXQ2  =  (A. 21) 

oqi  oq-2 

Note  the  right  hand  side  is  component  3  of  the  curl  of  the  electric  field  in  the  general 
coordinate  system. 

(V,  x  Ej  =  =  V  x  E  •  («!  x  u2)  Q,Q2  (A. 22) 

Similarly,  component  1  and  component  2  of  the  curl  of  the  electric  field  in  the  general 
coordinate  system  can  be  written  as 

/  \  1  fj  fp  fj  fp 

(v,x£j  =^-^-  =  \7xE-(u2xu3)Q2Q3  (A. 23) 

/  \  2  f)  TP  f)  fP 

( V,  x  E  Y  =  — A  -  —1  =  V  x  E  •  («3  x  u,)  Q3Q3  (A.24) 

V  /  oq3  oqi 

Note  how  the  curl  operation  in  the  general  coordinate  system  is  of  the  same  form  as  in 
Cartesian  coordinates,  with  the  only  difference  being  scale  factors  on  the  component 
parts  of  the  vector  field. 

The  form  for  the  curl  in  the  general  coordinated  system  can  be  substituted  into 
the  left-hand  side  of  Faraday’s  Law  (Equation  A. 22)  and  determine  Faraday’s  Law’s 
form  in  the  general  coordinate  system. 

/  ,  \  3  _  3H 

x  EJ  —  V  x  E  ■  (ui  x  u2)  QiQ2  =  •  («i  x  u2 )  gig2  (A.25) 
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The  magnetic  field  can  be  expressed  in  terms  of  its  contravariant  components  as 


H  =  H'in  +  H2u2  +  H3u3 


(A. 26) 


The  relationship  between  the  contravariant  and  covariant  components  of  a  vector  is 


1 

1 _ 

^  ui  ■  iii  iii  ■  u2  iii  ■  u3  ^ 
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u2  ■ iii  i 2  •  U2  «2  •  ^3 

to 

=  9 

H2 
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_ 1 

y  «3  •  Ul  M3  •  U2  U3-U3  J 

H3 

- 1 

CO 

(A. 27) 


where  g  is  the  metric  tensor  of  the  general  coordinate  system  (not  the  Euclidean 
metric).  A  more  compact  way  of  expressing  this  is 

3 

Hi  =  (A. 28) 

3= 1 

Solving  for  the  Hl  components  results  in 

3 

W  =  Hi,  (A. 29) 

3= 1 


where  g1]  are  the  components  of  g~l .  The  above  can  be  used  in  the  expressions 
to  solve  for  components  1,  2,  and  3  for  the  curl  of  the  electric  field  in  the  general 
coordinate  system. 

x  Ej  =  -/d0  ^  g^-^-iii  ■  (u2  x  u3)  Q2Q3  (A. 30) 

3= 1 

(vg  X  EJ  =  -Ido  ^2  g2j-g2-U2  ■  (u:i  X  u  1)  Q3Q1  (A.31) 

3= 1 

,  *  S  3  ^  QJJ 

\VqxEj  =  -/dp  g3j 3  •  (ui  x  u2)  Q1Q2  (A. 32) 

3= 1 
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By  defining 

/V1  =  fi3  \ui  ■  («2  x  u3)\Q1Q2Q3(QiQj)~1,  (A. 33) 

and 

Hj  =  QjH (A. 34) 
the  following  is  the  expression  for  the  components  of  the  curl  of  the  electric  field. 

(v9  x  Ej  (A. 35) 

3= 1 

A  similar  process  can  be  done  for  Ampere’s  Law  to  show  that 

(v,  *;/)•  =  ,„ (A. 36) 

3= 1 

Thus,  given  a  coordinate  transformation  from  Cartesian  coordinates,  the  behavior 
of  the  electromagnetic  fields  in  the  coordinate  transform  space  can  be  realized  in 
Cartesian  coordinates  using  a  complex  material  with  permittivity  and  permeability 
tensors  described  as  [97] 

eVJ  =  gVJ  |«i  •  («2  x  «3)|  QiQ2Q3{QiQj)-li  (A.37) 

jdlJ  =  gv  |mi  •  (u2  x  u3)\  QiQ2Q3(QiQj)~1-  (A.38) 
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Appendix  B.  Green’s  Function  Derivation 

A  Green’s  function  for  a  magnetic  line  source  radiating  in  the  presence  of  a  layered 
PEC  cylinder  is  derived.  The  geometry  for  this  problem  is  shown  in  Figure  B.l. 
The  final  solution  will  be  reached  by  first  solving  for  the  Green’s  function  for  a  PEC 


Observation  Point 


Magnetic  Line 
Source  in  Far  Zone 


Figure  B.l:  Problem  geometry  for  Green’s  function  derivation 

cylinder  illuminated  by  a  magnetic  line  source.  This  is  done  to  ensure  the  process 
is  correct  and  to  provide  a  series  of  checks  for  solution  accuracy.  Multiple  layers  of 
dielectric  materials  surrounding  the  PEC  will  then  be  added  to  arrive  at  the  final 
solution. 

For  this  problem,  the  source  is  an  infinite  magnetic  line  current  in  the  £  direc¬ 
tion.  Because  of  this,  the  incident  and  scattered  magnetic  fields  will  only  be  in  the 
£  direction.  Since  the  H  field  will  only  have  a  £  component,  the  vector  potential  F 
will  also  only  have  a  £  component;  the  vector  potential,  A  =  0.  The  potential  field 
must  obey 

V2F  +  k2F  =  -eM,  (B.l) 

which  is  a  vector  wave  equation.  Since  both  F  and  M  are  £  directed,  the  following 
scalar  equation  results 

V2Fz  +  k2Fz  =  -eMz.  (B.2) 
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Also,  z-invariance  has  been  assumed.  Thus, 


+  k2Fz  —  — eMz , 


(B-3) 


where  X72  is  the  transverse  Laplacian,  which  in  cylindrical  coordinates  is 


„2  1  d  (  d  \  Id2 

V2  =  (  r—  ]  + 


r  dr  V  dr  J  r2  dO2 


(B-4) 


The  boundary  conditions  are 


dFz 


dr 


=  0, 


(B.5) 


which  is  a  result  of  the  PEC  boundary  at  r  =  a.  Also  the  radiation  condition  must 
be  satisfied. 

I 

=  -jkFz  |r^oo  (B.6) 


dr 


Finally,  it  is  expected  that 


F{r,6)  =  F(r,9  +  2vr). 


(B.7) 


Now,  the  Green’s  function  must  solve 


V2G(r,r')  +  k2G(r,r')  =  —8(r  —  r'). 


(B.8) 


Note  that 


S(f  f0_4(r-r  >)6(e-ff) 


(B.9) 


and  that 


G(r,  f')  =  G(r,  9;  r',  9') 


(B.10) 


The  differential  equation  the  Green’s  function  must  satisfy  can  be  rewritten  as 


V2G(f,f')  +  k2G(f,f')  = 


5(r  -  r')S(9  -  9') 


(B.ll) 
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From  [25],  the  solution  can  be  found  by 


G(r,$;  r',0')  =  W  um(6)u’m(6')Gr(r, r':  A„,  I 


(B.  12) 


where  um  are  the  orthonormal  eigenfunctions  which  satisfy 


ggr,  T  —  0- 


(B-13) 


Also,  note  that 


k2  -  Xr  -  -i  =  0, 


(B.  14) 


and  that  Gr  satisfies 


ld_  (  dGA  ( ,2  _  M  G  =  S(r  -  r') 

r  dr  V  dr  )  \  r2  /  r  r 


(B.15) 


The  following  boundary  conditions  on  the  Green’s  function  are  enforced. 


dG(r,  9 ;  r',  9') 


(B.16) 


G(r,  9 ;  r',  0')  =  G(r,  9  +  2vr;  r',  0') 


<9G(r,  0;  r',  0') 
dr 


(B.  17) 
(B.18) 


By  having  G(r,9-,  r',9')  satisfy  the  same  boundary  conditions  as  Fz,  the  complemen¬ 
tary  solution  will  vanish. 


To  begin,  first  solve  for  um{9)  and  apply  the  appropriate  boundary  condition. 


Recall 


A  general  solution  for  um{9 )  is 


ggr,  T  Am,  Um(9)  0. 


(B.19) 


um(9)  =  A  cos  \f\m9  +  B  sin  V\m9. 


(B.20) 
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Applying  um(0)  =um{6  +  2n)  yields 


A  cos  VX m9  +  B  sin  \f\m9  =  A  cos  VAm(0  +  2n)  +  B  sin  \f\ m(9  +  2tt).  (B.21) 


The  only  way  this  equation  can  be  satisfied  is  if  V%n  =  m  where  m 
Therefore 

y/Xm  =  m,  m  =  0,  1,  2,  3..., 

A m  =  m2,  m  =  0,  1,  2,  3... 

The  general  solution  can  be  written  as 


0,1,2,.... 

(B.22) 


um{9)  =  A  cos  rnO  +  B  sin  m9,  m  =  0,  1,  2...  (B.23) 


A  cos  rri0  is  orthogonal  to  B  sin  m6.  Therefore,  it  is  known  that  A  cos  m6,  m  =  0, 1,  2... 
are  all  part  of  the  solution  of  eigenvectors.  Additionally,  Bsmm9,  m  =  1,2,...  are 
also  part  of  the  solution  for  the  eigenvectors.  Also  um{9)  has  been  defined  as  an  or¬ 
thonormal  set  of  eigenfunctions.  Therefore,  it  is  possible  to  calculate  the  coefficients, 
A  and  B  because 

27 r 

J  A2  cos2  m9d6  =  1,  (B.24) 

o 

2t r 

J  B2  sin2  m6d0  =  1.  (B.25) 

o 

Using  the  trigonometric  identities 


cos2  rriO  —  |  |  cos  2 m9, 
sin2  rnO  —  \  —  A  cos  2m9, 


(B.26) 


results  in  the  coefficients  being  found. 


2?r 

■  A2 

/  A2  cos2  m9d6  =  — 

27T  27T 

r  r 

d9+  cos  2  m9d9 

J  J 

J  2 

0 

.0  0 

(B.27) 
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This  simplifies  to 


However,  this  is  only  valid  if  m  ^  0.  Note  when  m 


(B.28) 

0,  the  integral  becomes 


2tt 


0 


1, 


which  yields 


Following  a  similar  process,  it  is  found  that 


B  = 


(B.29) 


(B.30) 


(B.31) 


There  is  no  issue  when  m  —  0  for  B  since  sin(0)  =  0  which  results  in  a  trivial  solution. 
To  make  things  easy  to  write,  express  A  and  B  as 


A,B=  — 
2vr 


(B.32) 


where 


Cm 


1  m  =  0 


2  m  =  1,  2,  3. 


(B.33) 


Note  this  shows  that  when  m  —  0,  B  —  y  This  isn’t  necessarily  true,  but  makes 
no  difference  since  the  term  is  multiplied  by  0  and  the  result  is  the  same.  Now, 
substituting  the  eigenvectors  results  in 


m= 0 


1 

27T 


OO 

em  [cos  rnd  cos  mQ'  +  sin  m6  sin  m6 '] . 

m= 0 


(B.34) 
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Using  the  trigonometric  identity 


cos  a;  cosy  +  sin  x  sin  y  =  cos  (a;  —  y), 


(B.35) 


the  result  is 

OO  1  OO 

E  Um(0)um(0')  =  —  E  em  cos  m(9  -  9'),  (B.36) 

m= 0  m= 0 

which  is  the  first  part  of  the  solution  for  the  Green’s  function. 

Next,  the  solution  for  G(r,  r';Xm)  is  found.  It  has  been  shown  \m  =  m2;  thus 
the  notation  can  be  rewritten  as  G(r,r';m2).  G(r,r';m2)  =  Gr  will  satisfy 


1  d  f  dGr\ 
r  dr  \  d  r  J 


(B.37) 


Multiplying  each  side  of  the  equation  by  r  results  in 


d  f  dGr  \ 

Tr  V  J  + 


— 8{r  —  r'). 


(B.38) 


This  equation  can  be  solved  using  the  U  —  T  method,  where  U(r)  and  T(r)  solve  the 
following: 


d_ 

dr 


dGr 

dr 


m 


+  k  r - 


u(r) 


=  0. 


(B.39) 


n  r) 


Note  U(r)  satisfies  the  boundary  condition  at  r  —  a  and  T(r)  satisfies  the  radiation 
condition  as  r  — >  oc.  The  above  equation  is  Bessel’s  equation,  and  Bessel’s  equation 
is  solved  by  Bessel,  Neumann,  and  Hankel  functions. 


U (r)  will  be  used  to  construct  the  solution  for  when  r  <  r'.  In  this  region,  either 
standing  waves  or  waves  radiating  outward  (depending  on  6)  are  expected.  Therefore, 
the  form  of  the  solution  for  U (r)  is  written  as 


U(r)  =  AJm(kr)  +  BH^(kr). 


(B.40) 
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The  boundary  condition  at  r  —  a  results  in  U'{r  =  a)  =  0.  Therefore,  the  relationship 
between  the  coefficients,  A  and  B,  can  be  found. 


B 


-A 


Uka ) 

H'£\ka) 


Letting  A  —  1,  the  result  is 


U  (r) 


Jm(kr) 


-40) 

Hm\ka) 


For  T(r),  the  form  of  the  solution  is 


T(r)  =  CH™(kr)  +  DH^(kr). 


T{r )  will  have  to  satisfy  the  radiation  condition.  This  requires  that 


T(r)  =  H^(kr), 


where  D  —  1.  The  Green’s  function  will  have  the  form 


Gr 


U(r<)T(r  >) 
c 


where  r<  is  the  lesser  of  r  and  r',  r>  is  the  greater  of  r  and  r',  and  c 
and  is  defined  as 


c  =  r[TU'  -  T'U}. 


The  relevant  equations  are 


u (r)  =  Jm(kr)  - 

U‘(r)  =  .VJkr)  -  ^ ^H2\kr), 
T(r)  =  HS\kr), 

T‘(r)  =  H‘2\kr), 


(B.41) 

(B.42) 

(B.43) 
7  =  0.  Thus 

(B.44) 

(B.45) 
is  the  conjunct 

(B.46) 

(B.47) 
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Solving  for  the  conjunct,  the  result  is 


c  —  r 


J'Jkr)H^(kr)  -  y 
Hm  (ka) 

-Jm(kr)H'S\kr)  HV\kr)H'M (At) 

Hrk'{ka) 


(B.48) 


which  simplifies  to 


c  =  r  [j'Jkr)H%\kr)  -  Jm(kr)H^\kr)\  .  (B.49) 


Knowing  that 

Hm(kr)  =  Jm(kr)  -  jYm(kr ), 
H'i2\kr)  =  JUkr)  -  fl”(kr), 

it  is  found  that 


(B.50) 


c  =  r  [J'm(kr)Jm{kr)^jJ'm{kr)Ym{kr)  -  J'Jkr)  Jm(kr)  +  j  Jm{kr)Y^{kr)} ,  (B.51) 
which  reduces  to 


c  =  jr[Jm(kr)Y^{kr)  -  J'm(kr)Ym(kr)}.  (B.52) 

Using  the  identity 

Jm(kr)Y'  (kr)  -  J'm(kr)Ym(kr)  =  — ,  (B.53) 

nr 

it  is  found  that 

c  =  j~,  (B.54) 

7T 


which  is  a  constant,  as  it  should  be. 

Next,  the  Green’s  function  is  written  as 

nSKkr>).  (B.55) 


r,  .7T 

Gr~~J2 


Jmikr <)  - 


Hm\ka) 


Hm  {kr 
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Subbing  back  into  the  standard  form 


G(r,9-  r',9')  =  ^  um(0)u*m(0')Gr(r,  r';  Am),  (B.56) 

m= 0 


the  result  is 


G(r,  6;  r' ,  O') 


-iE- 

m= 0 


,  cos  m{9  —  O') 


dm(fcr<) 


J'm{ka ) 
H'i2\ka) 


H{2\kr<) 


Hm\kr>). 


(B.57) 

To  ensure  this  is  the  correct  answer,  a  coordinate  shift  can  be  performed.  First,  let 
a  — »  0  and  then  shift  the  coordinate  system  such  that  r'  =  0  i.e.  move  the  origin  of 
the  system  to  the  location  of  the  line  source.  The  Green’s  function  should  reduce  to 
the  Green’s  function  for  the  radiation  from  a  magnetic  line  source  in  free  space. 

First  note  that  for  a  =  0,  the  term  Jji‘2{ka ^  will  vanish  for  all  m.  This  is  because 

Hm  ’(ka) 

for  all  m,  the  denominator,  Hm  (0)  goes  to  -oo.  Therefore,  if  a  — >  0,  the  Green’s 
function  becomes 


G(r,  0]  r',  O')  =  ^  emcos  m(9  -  Q')Jm(kr<)H^)(kr>).  (B.58) 

m= 0 


A  coordinate  transformation  can  performed  such  that  the  origin  is  shifted  to  r' . 


r<  =  0 


The  Green’s  function  can  be  rewritten  as 


(B.59) 


G(r,  9;  r',  O')  =  ^  emcos m{6  -  O')  Jm(0)H^)(k\r  -  r'\).  (B.60) 

m= 0 

It  is  known  that  Jm{ 0)  =  1  only  when  m  =  0.  For  all  other  possible  values  for  m, 
An(0)  =  0.  Thus,  the  summation  is  no  longer  required,  and  the  result  is 

G(r,6-y,ff)  =  -Iflf  (tlf-f'l),  (B.61) 
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which  is  the  correct  answer.  Thus,  the  methodology  used  is  correct,  and  the  Green’s 
function  for  an  infinite  magnetic  line  source  radiating  in  the  presence  of  a  PEC  with 
radius  a  surrounded  by  a  layered  material  can  be  developed.  For  the  geometry  shown 
in  Figure  B.l,  all  the  media  are  homogeneous.  Thus,  the  same  governing  wave  equa¬ 
tion  applies  with  the  only  change  being  to  the  wave  number.  Specifically,  the  equation 
that  must  be  solved  is 

Vf  Fz  +  k2Fz  =  —eMz  (B.62) 


where 


k  = 


k\  —  icUy/ /Arl^rll^o^-o  f  1  X  /  X  /' I 

^  k%  dlC J yj flri£rifl0£0  ri  <  V  X  V i-\-\ 

kn+l  —  k0  —  LO ^ J jJj0£o  T  T  n 


(B.63) 


where  i  —  2,  3,  ...n.  The  same  boundary  conditions  as  before  apply  and  are  repeated 


here  for  convenience 


=0 

dr  I  r=a  ’ 

fk|  =  _ jkF ; 

nr  r — x 


z  \r — xx)  5 


Fz(r,9)  =  Fz(r,0  +  2vr). 


(B.64) 


Additionally,  there  are  now  multiple  junction  conditions  at  r  =  r1;r  =  r2,  ...r  =  rn. 
These  are 

FI  =  F I  + 

z'r=rm  z\r=r+  J 

1  8FZ  I  _  1  diX  I 

£-  dr  \r=r^  e+  dr  \r=rti  ’ 

where  m  =  1,  2,  ...n  and  and  e+  are  the  relative  permittivities  in  the  regions  around 
the  boundary.  The  Green’s  function  solves 


(B.65) 


V2G(f,f')  +  k2G(f,f') 


5(r  -  r')S(9  -  9') 
r 


(B.66) 


The  solution  can  be  found  by 


G(r ,  9;  r',  9')  =  Um(0)u*m(9')Gr(r,  r';  Am). 

m= 0 


(B.67) 
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From  this  point,  the  development  proceeds  exactly  as  it  did  for  the  unlayered  PEC. 
It  has  been  shown  that 


^rn{0)u*m{6')  =  ^  €m  C°S  ~  9 l7)]  •  (B-68) 

m= 0  m= 0 


The  Green’s  function,  Gr(r,r Am),  can  be  found  using  the  U-T  method  where  U(r) 
and  T(r)  satisfy 


d 

dr 


dGt 

dr 


i,2  rn 

+  [  k2r - 

r 


U(r) 


=  0. 


(B.69) 


T(r) 


There  are  now  multiple  regions  necessitating  the  need  to  write  the  form  for  the  solu¬ 
tion  of  U{r)  and  T(r)  in  each  region  and  then  apply  the  junction  conditions  to  solve 
for  the  unknown  constants.  The  form  for  U (r)  in  each  region  is 


U(r) 


AlUhr)  +  ir)  a  <  r  <  ri, 

A%mJm{kir )  +  B'mH${kir)  i  <  r  <  rh 


(B.70) 


where  i  —  2,  3 ...n  +  1,  n  is  the  number  of  layers  material,  and  r  =  rn+\  is  understood 
such  that  rn+i  — oo.  Note  the  first  equation  must  satisfy  the  boundary  condition  at 
r  =  a.  B]m  has  already  been  found  in  previous  parts  of  this  exam,  and  it  is  the  same 
here  due  to  the  same  boundary  condition  at  r  =  a. 


■Uha) 

Hm\k\a) 


(B.71) 


A  similar  procedure  is  followed  when  solving  for  T(r).  First,  the  forms  of  the  solution 
for  T{r)  in  the  various  regions  are 


T(r) 


C'mJm(k,r)  +  D‘mH£Hktr) 


a  <  r  <  r1? 
n- 1  <  r  <  rt. 


(B.72) 


Like  before,  the  solution  for  the  coefficients  in  a  particular  region  are  already  known, 
except  in  this  case,  the  known  coefficients  are  C'^+1  and  D ”+1  due  to  the  same  bound- 
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ary  condition  as  r  — >  oo. 

C”+1  =  0 
B”+1  =  1 

Note  Gr  will  again  have  the  form 

^  f/(r<)r(r>) 


(B.73) 


(B.74) 


where  c  is  the  conjunct  and  is  equal  to 


c  = 


,;o  An+ 1 


IT 


(B.75) 


The  magnetic  line  source  will  always  be  in  the  free  space  region.  Additionally,  this 
formulation  will  only  be  used  such  that  r  <  r'.  Substituting  these  results  in  the  final 
form  for  the  Green’s  function,  G(r,  6 ;  r',  9')  for  a  PEC  cylinder  of  radius  a  surrounded 
by  n  layers  of  dielectrics  with  varying  radii,  r\  in  the  presence  of  an  infinite  magnetic 
line  source.  There  are  two  forms.  The  first  is  for  the  region  rj_i  <  r  <  r\. 


OO 

G(r,  9-,  r',  9')  cos  M*  ~  d')]  WMhr)  +  far)]  H%\k0r')  (B.76) 

4  m—0  Am 

The  second  is  for  the  region  r  >  rn. 

OO 

G(r,  9\  r',  O')  = -jT,  Th  cos  [m(9  -  9')\  [a”+1  Jm(fc0r<)  +  ^+1F^(fc0r<)  1  H^(kor>), 

4  m= 0  Arn 


where 


(B.77) 

A  n+ 1  _  1 

Jlm  -L> 

(B.78) 

(B.79) 
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and  i  =  2,3 ...n.  The  Alm  and  Blm  coefficients  are  found  by  applying  the  junction 
conditions  at  the  radial  interfaces. 


Al 


Al 


Uh r.)  -  3i*la)  gfffori) 

(kid) 


Hr, 


=  differ,)  +  BlH%\k2r,)  (B.80) 


-  AAr-»'S\k in) 

(Kia) 


=  ^(fcn)  +  BWN  (B.81) 

A2mJm{k2r2)  +  B2mH{2)  (k2r2)  =  A3mJm(k3r2)  +  BfnH{2)  (k3r2)  (B.82) 

AlJ'm{k2r2)  +  BlH'£\k2r2)  =  A3mJ'm(k3r2 )  +  B3mH'£\k3r2)  (B.83) 


AnmJm{knrn )  +  BnmH%\knrn )  =  Jm{k0rn )  +  Tr+1tf^(£;0rn)  (B.84) 

AnmJ'm{knrn)  +  BnmH'£\knrn)  =  J'm(k0rn )  +  ^+1^2)(fc0r„)  (B.85) 


These  equations  can  be  written  in  matrix  notation  of  the  form  Ar  =  B,  where 
A  is  a  2n  x  2n  matrix,  7?  are  the  forcing  functions,  and  x  is  the  solution  vector.  As 
an  example,  consider  a  case  where  four  layers  of  homogeneous  material  surround  the 
PEC  cylinder.  The  matrix  to  solve  for  the  unknown  coefficients  is 


H^\k0r3) 

H'ri2>(kor3) 

0 

0 

0 

0 


where  K, 

To  ensure  the  accuracy  of  the  derived  Green’s  function,  the  form  in  Equation 
B.77  was  used  to  determine  a2D  for  a  layered  cylinder.  The  results  were  compared 
to  those  obtained  using  a  Cornsol  simulation  for  a  simplified  cloak  with  material 


- Jr 

n(k3r3) 

-H^\k3r3) 

0 

0 

0 

"  BL  " 

~Jn 

%(,kar3) 

—  j'r 

n  (,k3r3) 

-H'S2\k3r3) 

0 

0 

0 

K, 

~Jn 

,(k0r3) 

Jrr 

i(k3r2 

H^\k3r2) 

—  Jm(k2r2) 

H^\k2r2) 

0 

Bm 

0 

J'h- 

t(k3r2 

(k3r2) 

-J!m(k2r2) 

( k2r2 ) 

0 

A2m 

0 

0 

0 

Jmifarx) 

Hrn)  (k2ri ) 

-Jm(feiri  +  KmH^  (kiri) 

Bm 

0 

0 

0 

Jm(k2ri) 

H&2)(k2r  i) 

—  J^ikiri  +  (fci^i) 

A\n 

0 

J'm(k  la) 
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FEM  and  GF  a2D  Comparison 


Figure  B.2:  Green’s  function  and  FEM  results  comparison 


parameters  put  forth  by  Yan  et  al.  and  shown  in  Equation  B.86. 


(B.86) 


In  order  for  the  Green’s  function  to  approximate  a  radially  varying  cloak  as  de¬ 
scribed  in  Equation  B.86,  the  number  of  layers  used  in  the  formulation  must  be  large. 
The  Green’s  function  results  were  determined  using  5,000  layers  to  approximate  the 
anisotropic  material.  The  FEM  results  were  obtained  with  MEL  =  0.01A.  The  calcu¬ 
lated  scattering  widths  from  the  two  methods  were  very  similar,  as  shown  in  Figure 
B.2.  The  A  for  these  results  is  0.004  m2,  which  is  quite  good.  There  is  a  noticeable 
difference  in  the  region  where  6  =  0°.  This  error  can  be  further  reduced  by  increasing 
the  number  of  layers.  Based  on  these  results,  it  was  concluded  the  Green’s  function 
was  correct. 
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Appendix  C.  Vector  Basis  Functions 

The  Comsol  Multiphysics  software  package  uses  the  finite  element  method  to  solve 
the  vector  wave  equation  and  associated  boundary  conditions.  The  formulation  used 
relies  heavily  on  that  developed  in  Jianming  Jin’s  book,  The  Finite  Element  Method 
in  Electromagnetics ,  Second  Edition  [33].  For  the  simulations  in  this  research  ef¬ 
fort,  tetrahedral  elements,  as  shown  in  Figure  C.l,  were  used  to  discretize  all  three- 
dimensional  domains.  Each  tetrahedral  element  has  four  vertices  with  coordinates 
(xf  ,yf,zf)  with  i  =  1,2,  3, 4  (black  numbers)  specifying  the  node  numbers.  Within 


each  tetrahedral  element,  the  electric  field  is  approximated  such  that 

E  =  E  FE]  (C-D 

3= 1 

_ k  e 

where  N  ■  are  the  vector  basis  functions  for  edge  j  (red  numbers)  and  the  Ee-  are  the 
unknown  coefficients.  Each  tetrahedral  element  will  have  four  nodes  and  six  edges. 

Vector  basis  functions  (or  edge  functions)  are  used  in  electromagnetics  because 
node-based  expansion  functions  are  not  able  to  accurately  represent  the  boundary 
conditions  associated  with  various  aspects  of  the  vector  fields  [95].  The  vector  basis 
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Table  C.l:  Edge  Node  Numbers 


Edge  i 

Node  ii 

Node  i2 

1 

1 

2 

2 

1 

3 

3 

1 

4 

4 

2 

3 

5 

4 

2 

6 

3 

4 

functions  effectively  enforce  the  requirement  for  continuity  of  tangential  fields  at  in¬ 
terfaces.  A  vector  basis  function  is  defined  for  each  edge  in  every  element.  Hence, 
each  element  has  six  vector  basis  functions.  These  vector  basis  functions,  vectN f,  can 
be  expressed  as 

N]  =  Wili2lt  =  [LhVLh  -  Lt2VLH\  If.  (C.2) 

Note  i  =  1,2,  ...6  and  is  the  edge  number  while  i\,2  =  1,2,  ...4  and  refer  to  node 
numbers.  An  edge  connects  node  i\  to  node  i2.  The  values  for  i\  and  i2  which 
correspond  to  the  edge  number  are  shown  in  Table  C.l.  As  an  example,  in  Figure 
C.l,  Ni  is  the  vector  basis  function  for  edge  1  (red),  which  connects  nodes  1  and  2 
(black).  Similarly,  Nq  is  the  vector  basis  function  for  edge  6  (red)  which  connects 
nodes  3  and  4.  The  If  term  is  the  length  of  edge  i. 

The  linear  interpolation  functions  (Lf,  i  =  1,2,  3, 4)  are  found  using  a  process 
developed  when  using  nodal-based  expansion  functions.  First,  the  linear  interpolation 
functions  can  be  written  as 


\x,y,z)  = 


6Ff 


( ai  +  bfx  +  cfy  +  dfz) . 


(C.3) 


The  yet-to-be-defined  terms  are  the  af,bf,cf,df  and  Ve  terms.  This  is  done  below. 

In  nodal-based  expansion  functions,  the  unknown  function  within  the  tetrahe¬ 
dral  element,  (pe ,  is  defined  as 


4>e(x,  y,  z)  =  ae  T  bex  +  cey  +  dez. 


(C.4) 
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Evaluating  0e  at  each  node  yields 


4>l  =  ae  +  bex  i  +  ceyi  +  dez\ 
0|  =  ae  +  bex2  +  cey2  +  dez2 
=  ae  +  feea;3  +  ce?/3  +  dez3 
04  =  ae  +  feeo;4  +  cey4  +  dez4 


(C.5) 


These  can  be  written  in  matrix  form  as 


^  1  x\  y{  z{  ^ 

1 

D 

3 

0? 

1  4  y2  ze2 

be 

02 

1  x%  y\  4 

ce 

0! 

v  1  x\  yl  zl  j 

de 

.  ^  . 

Cramer’s  rule  can  be  used  to  solve  for  ae,be,ce ,  and  de. 


(  4 

x\ 

?7i 

/i 

01 

y* 

4  ^ 

ae  =  ——  det 

01 

ry*  C 

1/1 

*1 

,  be=  )r  det 

1 

01 

i/l 

4 

6Ee 

01 

■4 

yl 

-3 

6Ee 

1 

0! 

yl 

4 

V  04 

ry*  C 

yl 

*4  ) 

V1 

04 

y% 

4  ) 

(l 

a;® 

01 

zl  ^ 

'l 

x\ 

i/i 

01  ^ 

e  1  i 

1 

/y>C 

x2 

01 

4 

Je  l  . 

1 

ry*  C 

yl 

0! 

c  =  — —  det 

,  d  =  — —  det 

6Ee 

1 

03 

4 

’  6Ee 

1 

®3 

yl 

03 

l1 

/y>  C 
*^4 

04 

4  ) 

V1 

/y>  C 

*^4 

yl 

04  ) 

Note  Ve  is  defined  as 


Ce 


/ 


V 


i 

i 

i 

i 


Xi 


nr*'-' 

0/0 


y\  4  ^ 
ye2  4 
y%  4 
yl  4  / 


(C.6) 


(C.7) 


(C.8) 


(C.9) 
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Table  C.2:  Example  Tetrahedral  Element  Node  Locations 


Node  i 

A 

vt 

A 

1 

2 

2 

2 

2 

1 

1 

1 

3 

3 

1 

1 

4 

2 

3 

1 

\Ve\  is  the  volume  of  the  element.  The  a?,  bp  Cp  and  dj  terms  are  found  by  evaluating 
the  determinants. 

ae  =  ^  +  ae2(pe2  +  a\<\)\  +  a|0|) , 

be  =  (6?0!  +  b\ (j)e2  +  b%<t>l  +  , 

1  (U.lUj 

ce  =  (cf^f  +  c|0|  +  c§0§  +  c|0|) , 

de  =  (dl<t>l  +  ^202  +  ^3^3  +  dl<l>l)  ■ 

Therefore,  the  aep  bp  Cp  and  d?  terms  are 


x2 

ye2 

*1 

A 

yf 

*i 

yf 

Ae 

2/i 

*i 

II 

4)  i-i 

<3 

x3 

2/1 

Z3 

,  «2  = 

x3 

2/1 

*3 

5  a3  — 

yf 

A 

,  «1  = 

2/1 

X% 

2/1 

zt 

x\ 

2/1 

*1 

X% 

2/1 

*1 

^3 

2/1 

*3 

1 

2/1 

Ze2 

1 

yf 

~e 
z 1 

1 

yf 

re 

21 

1 

yf 

zf 

61  = 

1 

2/1 

z3 

,  be2  = 

1 

2/1 

~e 

-3 

.  &3  = 

1 

2/1 

re 

z2 

,  &1  = 

1 

2/1 

z2 

1 

2/1 

z4 

1 

2/1 

z% 

1 

2 /! 

zt 

1 

2/1 

Z3 

1 

z2 

1 

*1 

1 

zf 

1 

7e 

Z1 

cf  = 

1 

x3 

Z3 

,  c2  = 

1 

^3 

Z3 

5  c3  = 

1 

Zf 

,  cl  = 

1 

*2 

7e 

Z2 

1 

X% 

z4 

1 

X% 

z4 

1 

£4 

zf 

1 

^3 

Z3 

1 

x\ 

2/1 

1 

xf 

yf 

1 

*! 

yf 

1 

xf 

yf 

d\  = 

1 

x3 

2/1 

,  <*£  = 

1 

x3 

2/1 

.  dt  = 

1 

2/1 

,  d!  = 

1 

x\ 

2/1 

1 

X% 

2/1 

1 

X% 

2/1 

1 

2/1 

1 

x3 

2/1 

In  the  above  matrices,  (xf ,yf,zf),  i  =  1,2, 3, 4,  are  the  coordinates  of  the  nodes  of 
the  tetrahedral  elements. 

As  an  example,  consider  a  tetrahedral  element  with  the  vertices  shown  in  Table 
C.2.  First,  evaluate  Ve. 
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(C.ll) 

Note  the  volume  of  the  tetrahedral  element  as  defined  in  Table  C.2  can  be  found  by 

v  =  lA°h  =  l(lw)1  =  l  <G12> 

Thus,  \Ve\  does  equal  the  volume  of  the  element.  Next,  evaluate  the  a®,  6®,  c®,  and 
d®  elements  using  the  matrices  defined  above. 


a\  =  4,  a  2  —  —8,  a \  =  0,  aeA  =  0, 
b\  =  0,  be2  =  2,  61  = -2,  61  =  0, 
Ci  =  0,  =  1,  c|  =  1,  =  -2, 

d^  =  -4,  d2  =  1,  dg  =  1,  d\  =  2. 

Using  these  values,  the  linear  interpolation  functions,  L®,  are 


(C.13) 


L{(x,y,z)  =  -l  +  z, 

Le2(x,  y,z)  =  ~l(2x  +  y  +  z-8), 
Le3(x,y,z)  =  -\(-2x  +  y  +  z), 
L\{x,y,z )  =  -\{-2y  +  2z) , 

with  the  gradients  being 


(C.14) 


VL®  =  5, 

VT®  =  -|(2x  +  £  +  z), 
VLs  =  -i(_2£  +  £  +  *)* 

VLt  =  -\(~2y  +  2z). 


(C.15) 
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The  lengths  of  each  side  in  the  tetrahedral  element  are 


It  =  V3,  Z|  =  n/3,  ^3  =  \/2, 

ZI  =  2,  Zf  =  V5,  /6e  =  a/5- 

It  is  now  possible  to  determine  the  N?  functions. 

W12  =  -\  [(2z  -  2)  x  +  (z  -  1)  y  +  (7  -  2x  -y)  z\ , 

^13  =  [(2  -  2z)  x  +  (*  -  1)  y  +  (2x  -  y  -  1)  z] , 

^i4  =  -i[(2-2*)y+(2j/-2)  i], 

^23  =  [{y  +  Z  +  4)  x  +  (~x  +  2  )y  +  (■ -x  +  2)z], 

W42  =  -\[{y—z)x  +  (-x  -  z  +  4)  y  +  (x  +  y  -  16)  z] , 
^34  =  -\ [(y  ~  z)  x  +  (-x  +  z)  y  +  (x  —  y)  z] . 


(C.16) 


(C.17) 


Note  that  for  all  W.^,  the  following  hold: 


V  •  Wili2 
v  x  Wili2 
ei  •  VL?  = 
e*  •  VL?  = 


0, 

=  2VLfx  x  VL?a, 


(C.18) 


where  e_,-  is  a  unit  vector  for  edge  i  pointing  from  node  i\  to  node  i-i-  By  satisfying 
these  constraints,  all  Wlxl2  vectors  have  a  constant  tangential  component  along  the 
ith  edge  while  having  no  tangential  component  along  the  five  other  edges.  Thus,  the 
vector-based  edge  elements  are  well-suited  for  electromagnetics  problems.  The  final 
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vector-based  edge  elements  are 


N,  =  W12lt  =  [(2 z  -  2)  x  +  (z  -  1)  y  +  (7  -  2x  -  y)  z] , 

N2  =  W13l%  =  [(2  -2z)x  +  (z-l)y  +  (2x-y-l)z], 

N3  =  Wule3  =  -^  [(2  -  2 z)y  +  (2 y  -  2)  2] , 

iv4  =  =  -\ [(y  +  z  +  4)  X  +  (-X  +  2)  y  +  (-X  +  2)  i] , 

iV5  =  kh42/§  =  [(y  -  z)  x  +  (-x  -  z  +  4)y  +  (x  +  y  -16)z], 

N6  =  Wul*  =  [(y  -  z)  x  +  (-x  +  z)  y  +  (x  -  y)  z] . 


Thus,  given  any  tetrahedral  element  within  a  domain,  the  above  method  is  used  to 
determine  the  linear,  vector-based  edge  elements  which  will  approximate  the  solution 
within  each  element.  This  is  done  for  each  element  within  the  domain.  The  unknowns 
are  manipulated  into  a  system  of  equations  which  are  solved  using  standard  techniques 
such  as  conjugate  gradient. 
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